1/2 


MICSOCOP''  RESOLUTION  test  CHART 

NA*  CS*L  - 'B65  '- 


AD- A 158  973 


^  Uuclassified 


M  security  classification  of  this  page 


1«.  REPORT  SECURITY  CLASSIFICATION 

Jnclasslfled 


L«-«ECURITY  CLASSIFICATION  AUTHORITY 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


3b.  OECLASSIFICATION/OOWNGRAOING  SCHEDULE 


*.  PERFORMING  ORGANIZATION  REPORT  NUM8ERIS) 


6a  name  of  PERFORMING  ORGANIZATION 


b.  OFFICE  SYMBOL 
(IfapplieabUt 


University  of  Pittsburgh 


6c.  ADDRESS  (Cily.  Stale  and  Z/P  Code) 

Department  of  Electrical  Engineering 
University  of  Pittsburgh 
Pittsburgh,  PA  15261 


sa  name  of  funoing/sponsoring 

ORGANIZATION  Air  Force  Office 

of  Scientific  Research 


8b.  OFFICE  SYMBOL 
Ilf  appUeable) 


.AFOSR/NM 


Sc.  address  iCity.  Stale  and  ZIP  Code) 

Division  of  Mathematical  and  Information 
Sciences;  Bolling  Air  Force  Base 
Washington,  DC  20332 


11.  title  (Include  Security  Clasttficationl 

hift-Variant  Multidimensional  Systems 


13.  personal  AUTHOR(S) 

^.'rofessor  N.K.  Bose 


TYPE  OF  REPORT 

Final 


16.  SUPPLEMENTARY  NOTATION 


3.  DISTRIBUTION/AVAILABILITY  OF  REPORT 

Approved  for  Public  Release  - 
Distribution  Unlimited 


S.  MONITORING  ORGANIZATION  REPORT  NUMBERISI 

AFOSR  TR-  7  ?  4 


7a  name  of  MONITORING  ORGANIZATION 

Air  Force  Office  of  Scientific  Research 


7b.  ADDRESS  ICity.  State  and  ZIP  Code) 

Division  of  Mathematical  and  Information 
Sciences;  Bolling  Air  Force  Base 
Washington,  DC  20332 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

Department  of  the  Air  Force 
Grant  No.  AFOSR-83-0038 


10.  SOURCE  OF  FUNDING  NOS. 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

ELEJWENT  NO.  . 

NO. 

NO. 

NO. 

2304 

A6 

13b.  TIME  COVERED 

FROM  2/1/83  TO  3/31/8 


14.  DATE  OF  REPORT  lYr..  Mo.,  Day) 

5/29/85 


15.  PAGE  COUNT 


COSATI  COOES 


GROUP  SUB.  GR. 


18.  SUBJECT  TERMS  iContInue  on  reverte  if  neceieary  and  identify  by  block  number) 

Multidimensional  systems  State-space  model 
Spatio-temporal  processing  Image  Restoration 
Shift-variant  systems  Bilinear  degradation 


19.  ABSTRACT  iConllnue  on  reverte  if  neeeitory  and  Identify  by  block  number) 

Procedures  for  compensating  for  effects  which  degrade  the  accuracy  of  remotely 
sensed  datas  by  mathematically  inverting  some  of  the  degrading  phenomenon  are  required  in 
biomedical,  industrial,  surveilance  and  earth  and  space  applications.  Images  to  be  re¬ 
stored  are  often  degraded  by  a  linear  spatially  varying  operation.  Degradations  due  to 
motion  blurring  and  optical  system  distortions  often  require  the  imaging  systems  to  be 
analyzed  by  modeling  the  degradation  as  a  linear  shift-invariant  motion  blur,  but  coun¬ 
terparts  of  such  techniques  in  the  shift-variant  case  is  severely  handicapped  by  the 
required  space-time  computational  complexity.  Furthermore,  a  systematic  approach  for 
handling  the  analysis  of  multidimensional  shift-variant  systems  remains  to  be  developed  in 
view  of  some  of  the  basic  differences  between  the  properties  of  shift- invariant  and  shift- 
variant  systems.  For  example,  the  crucial  factor  of  determining  the  best  sampling  rate 
in  digital  implementations  of  analog  shift-variant  systems  is  considerably  complicated  by 
the  fact  that  for  the  shift-variant  case,  the  output  bandwidth  can  exceed  the  input  band- 


31.  ABSTRACT  security  CLASSIFICATION 


Unclassified 


33b.  TELEPHONE  NUMBER 
Uylude^rnaPoiieJ^- 


DD  FORM  1473,  83  APR 


EDITION  OP  1  JAN  73  IS  OBSOLETE. 


SfiCURITY  CLASSIPICATION  OF  THIS  PAGE 


Unc  i  I  1  eii 


SECURITY  CLASSIFICATION  CF  THIS  FACE 


width. 

The  research  being  reported  was  directed  to  the  alleviation  of  the  shortcomings 
referred  to  above.  The  following  specific  results  have  been  obtained.  For  any  2-D  discrete 
first  quadrant  quarter-plane  causal  linear  shift-variant  (LSV)  system,  whose  Impulse  res¬ 
ponse  is  a  K-th  order  degenerate  sequence,  a  K-th  order  state-space  model  was  obtained. 

This  model  is  recursive  and  is  based  on  a  three- term  recurrence  formula  relating  any  point 
in  the  state-space  to  its  three  closest  neighboring  points  and  the  current  input.  The 
state-space  model  was  extended  in  order  to  model  2-D  discrete  LSV  systems  with  support  on  a 
causality  cone.  Subsequently,  the  2-D  quarter-plane  causal  and  weakly  causal  discrete 
models  were  generalized  to  the  n-D(n>2)  case.  The  resulting  state-space  models  are  recur¬ 
sive  and  are  based  on  a  (2^-l)-pointa  recurrence  formula,  which  for  the  causal  case  uses  the 
(2"-l) -closest  neighboring  points  in  addition  to  the  input  in  order  to  compute  any  current 
output  state.  For  the  weakly  causal  case,  the  (2''-l)  computed  outputs  required  are  not,  in 
general,  the  closest  neighbors  to  the  output  presently  being  computed.  Conditions  for  the 
existence  of  a  2-D  state-space  model  for  the  Inverse  system  were  obtained  and  with  these 
conditions  satisfied,  a  state-space  model  of  Che  inverse  system  is  readily  derivable  from 
the  original  one.  Models  for  the  2-D  LSV  system  and  its  inverse  can  be  used  to  perform 
analysis  and  deconvolution  problems  very  efficiently.  This  was  substantiated  from  derived 
expressions  for  space-time  computational  complexities. 

* 

Examples  of  physically  motivated  applications  making  use  of  the  theoretical  results 
developed  have  been  worked  out.  These  applications  include  effects  of  1-D  LSV  motion  blur 
and  the  blurring  due  to  Seidel  aberration  of  a  lens;  in  particular,  the  2-D  LSV  coma  aber¬ 
ration  was  studied  in  detail.  The  reconstruction  of  the  original  object  from  the  LSV 
blurred  image  was  carried  out  successfully  by  means  of  the  state-space  model  for  the  in^ 
verse  system.  For  the  construction  of  the  state-space  model,  the  impulse  responses  of  t».»r 
blurring  phenomena  were  approximated  in  a  degenerate  form  via  series  expansion  using 
orthogonal  functions. _  _ 


The  problem  of  restoring  Images  degraded  by  phenomena,  which  can  be  modeled  by  linear 
shift-invariant  systems,  has  been  discussed  above.  Sometimes,  the  degrading  phenomena  may 
not  be  accurately  modeled  by  systems  that  are  restricted  to  be  linear.  In  such  cases,  the 
incorporation  of  the  second-order  term  of  a  Volterra  series  (characterizing  the  input/output 
behavior  of  the  nonlinear  system),  which  forms  a  class  of  bilinear  systems,  improves  sub¬ 
stantially  the  measure  of  adequacy  for  the  model.  This  type  of  bilinear  syatem  occurs  iit 
imaging  through  turbulent  atmosphere,  coherent  Imaging  through  systems  with  time-varying 
pupils,  etc.  Recent  research  results  on  bilinearly  distorted  images  are  available  based  on 
finite  impulse  response  linear  digital  filtering  and  Bayesian  methods.  Here,  additional 
methods  have  been  investigated  for  restoring  images,  wliich  are  distorted  by  a  system  that  is 
describable  by  the  second-order  term  of  the  Volterra  series.  When  the  blurring  phenomenon  is 
nonlinear  and  can  be  modeled  either  by  a  shift-invariant  of  shift-variant  bilinear  system, 
the  data  restoration  problem  can  be  most  conveniently  formulated  as  a  special  system  of 
linear  equations  with  nonnegative  coefficients  whose  solution  is  required  to  satisfy  con¬ 
straints  like  nonnegativity  in  addition  to  being  factorable  with  the  factors  having  a  certain 
characterizing  property.  An  algorithm  implementing  this  objective  along  with  another 
important  alternative  applicable  to  a  specialized  model  are  discussed  in  the  research  report. 
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To  a  great  extent  the  techniques  for  analysis  and  res¬ 
toration  of  images  has  been  developed  under  the  assumption 
that  the  system  is  linear  shift-invariant  (LSI).  These 
techniques  are  successful  in  some  cases  because  a  system 
which  is  diffraction-limited  or  a  system  whose  object  plane 
undergoes  uniform  linear  motion  perpendicular  to  the  system 
reference  axis  does  indeed  satisfy  these  assumptions. 

However,  LSI  systems  are  singled  out  for  study  mainly  be¬ 
cause  of  the  widespread  understanding  of  the  Fourier  Trans¬ 
form  theory  along  with  well-known  fast  algorithms  for  its 
implementation.  In  comparison  with  LSI  systems,  very  little 
work  has  been  done  on  linear  shift-variant  (LSV)  sy s terns  .  Most 
of  the  research  on  two  dimensional  (2-D)  LSV  systems  has  been 
done  on  restoration  techniques  by  means  of  coordinate  trans¬ 
formations.  This  technique,  decomposes  the  LSV  system  into  a 
distortion  of  the  input  plane  followed  by  a  s h i f t- inva riant 
operation  and  terminated  by  a  distortion  of  the  output  plane. 
Essentially,  the  shi f t- va r ian t  problem  is  transformed  into  a 
sh i f t - inva r ian t  one  and  the  deconvolution  or  inverse  filtering 
is  done  with  the  well-known  sh if t - inva r i an t  methods  for  this 
purpose.  The  main  drawback  of  image  restoration  by  coordinate 
transformation  is  that  it  can  only  be  applied  to  a  limited  class 
of  LSV  systems.  For  one  dimensional  (l-D)  LSV  systems,  there 
has  been  some  research  activity  and  most  of  the  work  has  to  do 
with  the  reconstruction  of  motion  degraded  images.  Besides  the 
already  mentioned  approaches  to  image  restoration  based  on 


coordinate  transformations,  some  methods  have  used  1-D  state-space  models  for 
a  class  of  LSV  systems.  At  present  the  technique  for  analysis  and  deconvolu¬ 
tion  of  LSV  systems  are  either  applicable  to  the  1-D  case  or  to  a  very  restric¬ 
tive  class  of  2-D  LSV  systems.  It  is  important  to  note  that  Xhe  characteriza¬ 
tion  of  a  LSV  system  by  means  of  the  superposition  integral  or  the  superposi¬ 
tion  sum  holds  in  general,  but  it  is  completely  impractical  for  the  analysis 
and  the  deconvolution  of  LSV  systems.  Therefore,  we  can  conclude  that  there 
is  a  great  need  for  an  efficient  and  convenient  model  for  LSV  systems  which 
will  permit  the  analysis  and  deconvolution  in  a  simple  form. 

One  of  the  primary  objectives  of  this  research  has  been  to  provide  not 
only  a  mathematical  structure  for  the  state-space  modeling  of  discrete  LSV  sys¬ 
tems  but  to  apply  this  model  to  the  problems  of  efficient  analysis  and  decon¬ 
volution  of  multidimensional  systems.  It  is  expected  that  the  state-space 
model  will  be  useful  in  system  analysis  and  synthesis.  The  state-space  model 
should  be  useful  in  solving  the  problem  of  deconvolution  very  efficiently.  In 
fact,  for  the  noise-free  case  the  deconvolution  scheme  developed  from  the  state- 
space  model  is  considerably  more  efficient  from  space  as  well  as  time  computa¬ 
tional  complexity  standpoints  than  other  techniques  which  are  currently  avail¬ 
able. 

A  portion  of  the  time  for  the  proposed  research  was  also  devoted  to  the 
approximation  problem  of  specified  2-D  linear  shift-variant  impulse  responses 
by  K-th  order  degenerate  approximants . 

Several  applications  require  the  restoration  of  bilinearly  degraded  images. 
For  example,  in  coherent  optical  image  processing,  the  bilinear  term  in  the 
Vol terra  series  expansion  plays  an  important  role  due  to  the  nonlinear  nature 
of  partially  coherent  image  formation.  The  problem  of  restoration  of  images 
distorted  by  nonlinear  nonzero-spread  systems,  which  appear  quite  frequently 
in  real  world  problems  have  not  been  systematically  tackled.  Therefore,  portion 
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of  the  research  effort  was  directed  to  the  development  of  algorithms  for  digital 
restoration  of  bilinearly  degraded  images  (shift-invariant  as  well  as  shift- 
varian^ 
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c.  Details  of  Research  Results  Obtained 


The  primairy  outcomes  of  the  research  effort  are  in  the  following 
two  categories: 

1.  Linear  shift-variant  multidimensional  systems. 

2.  Restoration  of  bilinearly  degraded  images. 

Detailed  documentation  of  results  obtained  in  the  above  mentioned  areas 
are  contained  in  the  succeeding  pages  of  this  report. 
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LINEAR  SHIFT-VARIANT  ^^ULTIDIMENSIONAL  SYSTEMS 
5 . 1  Introduction 

The  study  of  1-D  time-varying  systems,  now  we 11- documented 
in  the  literature,  was  motivated  to  a  great  extent  by  the 
need  to  design  adaptive  control  systems,  where  the  parameters 
of  the  controller  are  adjusted  to  counterbalance  the  fluctu¬ 
ations  in  process  dynamics  stemming  from  changes  in  environ¬ 
ment.  For  example,  variations  of  flight  conditions  (including 
flight  speed)  of  supersonic  aircrafts  and  missiles  during  rapid 
ascent  through  the  atmosphere  introduce  time-varying  parameter 
variations.  In  space  flight,  relevant  problems  that  had  to 
be  tackled  included  those  [5.1]  of  transferring  a  space  vehicle 
from  one  orbit  to  another,  rendezvous  and  interplanetary 
guidance.  The  development  of  the  rigid-body  equations  of 
motion  of  artificial  satellites  requires  accounting  for  cyclic 
variations  of  inertial  torques  that  occur  as  the  vehicle 
progresses  in  its  orbit.  Similarly,  successful  solution  of  the 
satellite  rendezvous  problem  requires  dealing  with  the  combined 
complexity  of  both  mass  and  orbital  variations  ,  for  with 
both  the  target  and  interceptor  in  orbit,  the  interceptor 
must  expend  fuel-mass  in  maneuvering;  further,  the  kinematics 
of  interceptor  guidance  introduce  additional  time-varying 
parameters  forcing  the  characterizing  differential  equations 
to  have  both  periodic  and  aperiodic  coefficients.  Also  or 
interest  is  a  class  of  problems  which  relate  to  the  gyroscopic 
stabilization  of  orbiting  satellites.  In  circuit  and  systems 
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theory  specific  applications  concerned  with  the  theory  of 
modulators,  parametric  amplifiers  and  harmonic  generators  use 
linear  circuits  with  periodically  varying  parameters  [5.2]. 

The  developed  tools  of  multidimensional  systems  theory, 
as  documented  in  [5.3],  provide  scopes  for  considerable  appli¬ 
cations  in  a  broad  range  of  problems.  However,  over  the  last 
few  years,  it  is  fair  to  say  that  the  theory  of  linear 
shift-invariant  (LSI)  two  or  more  dimensional  systems  have 
been  well  understood.  However,  with  the  increased  activities 
in  the  area  of  image  processing  (optical  and  digital)  many 
physical  problems  are  charaterizable  by  a  2-D  linear  integral 
operator, 

»  oo 

g(x,y)  =  ;  /  f(u,v)  h(x,y;u, v)dv 

—  00  ••  oo 

where  g(x,y),  f(x,y)  are,  respectively,  the  output  and  input 
while  h(x,y;u,v)  is  the  system  response  at  (x,y)  to  an  unit 
impulse  applied  at  (u,  v)  .  The  presence  of  additive  noise  in 
the  preceding  input-output  model  of  a  physical  process  could 
also  be  expected.  In  1977,  Goodman  state  that  "the  type  of 
coherent  optical  processing  which  is  by  far  the  least  explored 
to  date  is  that  of  linear  space-variant  filtering"  [5.4], 
and  the  importance  of  overcoming  this  deficiency  was  realized 
because  space- variant  processing  is  required  in  various 
optical  and  digital  data  processing  applications  including 
restoration  of  images  degraded  either  by  space-variant 
aberrations  or  space- variant  motion  blur,  restoration  of  radio¬ 
graphs  blurred  by  a  space- variant  source  penumbra,  restoration 
of  rotation  blur  of  2-D  patterns  where  different  objects 
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suffer  from  different  blur  according  to  their  distance  from 
the  center  of  rotation  and  aldo  their  angular  positions  [5.5], 
and  performance  of  particular  transformations  (like  the  Mellin 
transform  and  Abel  transform) .  Also,  emerging  applications  in 
areas  like  robotic  vision,  have  necessitated  the  development 
of  techniques  suitable  for  recognition  patterns  which  may  be 
subjected  to  the  effects  of  unknown  rotation  and  scaling.  One 
such  technique  is  based  on  the  use  of  scaled  transform  (a 
modified  Fourier  transform) ,  which  have  been  shown  to  be 
equivalent  to  a  class  of  linear  shift- variant  (LSV)  filters, 
referred  to  as  the  class  of  "form- invariant"  filters  [5.6]. 

Actually,  the  problem  of  2-D  images  degraded  by  spatially 
varying  point  spread  functions  has  been  tackled  since  about 
the  early  seventies.  Sawchuk  [5.7]  transformed  the  spatially 
varying  problem  to  a  spatially  invariant  one  via  coordinate 
transformations  and  as  a  result  some  generality  in  the  approach 
had  to  be  sacrificed,  especially  with  the  presence  of  additive 
noise.  In  that  situation,  however,  statistics  of  the  noise  and 
object  random  processes  could  change  under  a  geometrical 
coordinate  transformation.  The  usual  assumption  of  the  rainimuir. 
mean- square  error  estimation  (MMSE)  method  is  that  these 
statistics  are  jointly  stationary,  and  this  may  not  be  valid 
in  many  cases.  Frieden  [5.8]  developed  a  positive  (an  incoherent 
object  scene  is  a  spatial  radiance  distribution  and  cannot 
be  negative)  restoring  formula  from  a  statistical  communication 
theory  model  of  image  formation  by  which  the  optimal  restored 
object  from  a  set  of  candidate  objects  is  required  to  obey  a 
principle  of  maximum  entropy.  Frieden 's  approach,  though 
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The  particular  structure  of  the  matrices  in  (5.34)  to  (5.36) 
will  allow  us  to  perform  the  summation  in  (5.32)  by  means  of 
a  (2  -1) -point  recursion.  This  recursion  will  be  a  gener¬ 
alization  to  k-dimensions  (k  >  2)  of  the  recursion  in  (5.10) 
and  the  justification  will  be  developed  next.  To  do  this, 
the  following  definitions  are,  first,  introduced. 

Definition  5.3:  The  .]Li  (n.+l)  row- vector  6,.  is 

J-i  j  — K, 

defined  as 


^1^2' 


A  [0  0 


I2I3 

0^-1 


0  .  .  •  0  ]  .  (5.37) 


(i,+l)-th  block  position 
from  the  right 


where  0  is  a  tt  (n.+l)  row-vector  of  zeros,  there  are  n, 

j=2  J 

such  0  vectors,  and  the  position  of  the  remaining  nonzero 
row-vector  is  indicated  above.  Furthermore,  the  (n^+1) 

K 

_ _ _ _ _ r  _ Z _ 


row-vector 


is  defined  as 


4  [0  0 


0  1  0  . . .  Oj 

I 


(5.38) 


(i, +l)-th  position  from 
the  right 


Definition  5  .  A  :  ^(j)  L  sum  of  all  (.)  distinct  vectors 

,^1^2-  •  -^k  ^  ^  .0 

^  i  '  or  ( 1=1 , 2  ,  .  .  .  ,k)  and  i^+i2+.  •  .+i^=j 


definition  5.5:  4  distinct  vectors 
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where  .  ,n^)  .  .  .  ,K)  and  s^(m^ ,  .  .  .m^^)  (i=l,2 . 

K)  are,  respectively,  linearly  independent  functions  of 
(n^  ,  .  .  .  and  (in^  iti^  •  •  •  •  ■ 

Substituting  (5.30)  in  (5.29).  after  carrying  out  some 
manipulations  analogous  to  the  ones  in  Section  5.2,  we  have 


y(n^,  .  .  .  ,nj^)  =  z^a^(n^,  .  .  .  ,n^)x^(n^,  .  .  .  ,n^),  (5.31) 


where 
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Xj^(n^,  .  .  ,n^)  =  E  ...  Z  8^(m^,  .  .  ,mj^)  •u(m^,  .  .  ,mj^).  (5 . 32) 

m.,  ®*0 

^  **  i=l,2,...,K 

The  expression  in  (5.32)  can  be  written  in  matrix  form 

3^  *  ^  ^^1  *  ’  *  ’  i“l  »2,..,K,  (5.33) 

where  x^  and  ^  are  vectors  of  order  ( [  (ti^^+l)  •  (^2+1)  •  •  • 

(n^+1) ]  X  1)  obtained,  respectively,  by  ordering  the  k- 

dimensional  arrays  {x^(i^,  i2 . i^)}  and{u(i2^,  i2.  •  •  • 

(i^“0, 1,  .  .  .  ,nj  ;  j=l,2,  .  .  .k)  by  means  of  a  lexicographic  ordering, 
described  next.  The  element  u(j^, j2 . • • • » occurs  in  row 
k-1  k-1 

i=l  ^  ii=i  ^ 


for  j^=»0 , 1 ,  .  .  .  ,n^ ,  i=l,2,...,k,  of  the  vector  A 

similar  lexicographical  ordering  is  given  to  the  elements 
tx^(j2_.  32’ '  •  •  ’  ^k^  '  vector  i=l,2,....K. 

-i  k^^l’ ■ ■ ’^k^  (5.33)  is  a  square  matrix  of  order 

( (n^+1)  •  (n2+l)  .  .  .  (t\+1)  )  .  given  by 


5.13 


into  its  components  and  summing  over  them  in  equation  (5.9). 
This  derivation,  even  though  simple,  cannot  be  easily 
extended  to  higher  dimensions.  In  this  section,  the 
particular  structure  of  a  matrix  will  be  used  to  our  advan¬ 
tage  in  obtaining  a  k-D  recursive  model.  First,  this  will 

be  done  for  the  case  when  the  impulse  response  sequence 
has  support  in  a  positive  cone  and,  subsequently,  weakly 
causal  LSV  systems  will  be  tackled. 

5.3.1  k-D  |*ositive  Cone  Causal  State-Space  Model 

A  k-D  positive  cone  causal  discrete  LSV  system  can  be 
characterized  by  the  superposition  svnn 


’"l 

2.  ...  Z  hCn^ ,  .  .  , n^^ ; m^ ,  .  .  ,m^)u(m.^,  .  .  ,mj^)  , 

m,  =*0  m.  =0 

^  (5.29) 


where  u(m)  is  the  input  at  coordinate  point  m=(m^  ,m2 ,  .  .  .  • 

y(n)  is  the  output  at  coordinate  point  n  =  (n^,n2 , . . . ,n^) 
and  h(n,m)  is  the  response  of  the  system  at  the  point  n 
to  a  unit  impulse  at  the  point  m.  At  this  point,  a 
definition,  which  is  a  natural  extension  of  Definition  3.1, 
is  given. 

Definition  5.2:  A  sequence  {h(n^,n2,  .  .  ,n^:m^,m2,  .  .  ,mj^) } ,  is 
K-th  order  degenerate  if  it  can  be  expressed  as 

K 


•■“k) 


r  aj^(ni,..,ti^)3i(!ni - (5.30) 

1=1 


h(n^, . . .n^;m^ 


I  • 
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x(n^,n2)  =  Aj^(n^,n2)x(n^-l,n2)  +  A2(n^.n2)x(n^,n2-1) 

(5.22) 

-  A2(n^,n2)x(n^-1, 112-1)  +  b(n^,n2)u(n^,n2)  , 

y<n^,n2)  =  c(n^,n2)x(n^,n2) ,  (5.23) 

where 


A^^  (n^ ,  n2) 

—  T(n^,n2)  ‘  T  (n^-l,n2) , 

(5.24) 

A2  (^2_ » ^2  ^ 

*  T(n2^,n2)  ■  T  (n^,n2-l)  , 

(5.25) 

A2(n^,n2) 

=  T(nj^,n2)  •  T*^(nj_-l,n2-l)  , 

(5.26) 

A 

f  112) 

=■  T(n^.n2)b(n^,n2)  , 

(5.27) 

£(n^,n2) 

»  c(nj_,n2)T”^(n^,n2)  . 

■(5.28) 

From  the  above,  we  can  say  that  the  state-space 
representation  in  (5.22)  to  (5.28)  is  related  to  the  state- 
space  model  in  (5.12)  to  (5.16)  by  the  transformation 
matrix  in  (5 . 21)  . 

5.3  k-D  State-Space  Model 

Now,  the  results  obtained  in  sections  5.2  will 
be  extended  to  the  k-dimensional  (k  >  2)  case.  The  state- 
space  model  in  Section  5.2  was  obtained  by  making  use  of  a 
geometric  argument  which  consists  of  splitting  a  2-D  grid 
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x(-l,n2)  =  0  for  tl2  =  0,1,..., ^2’  (5.20) 

The  state-space  model  introduced  in  (5 . 12)- (5 . 16)  is  based 
on  a  three  term  recursion,  which  means  that  in  order  to  gen¬ 
erate  the  state  vector  x(nj^,n2)  at  a  point  (n2^,n2)  we  need 
to  have  previously  computed  the  state  vectors  at  the  points 
(nj^-l,n2)  ,  (n^,n2-l)  and  (n^-l,n2-l)  in  addition  to  the 
current  input.  See  Figure  5.2. 

From  (5.10),  it  is  easy  to  see  that  the  states  are 
decoupled;  this  is  a  nice  feature  which  permits  parallel 
computation  of  the  states. 

The  state-space  model  introduced  here  has  a  structure 
similar  to  the  Fomas ini-Marches ini  state- space  representa¬ 
tion  15.3].  The  main  difference  between  these  two 

models,  is  that  in  (5,12)  and  (5,13)  some  of  the  matrices 
are  non- constant. 

The  state-space  model  in  (5.12)  to  (5.16)  is  not  a 
unique  state-space  model  representation  of  the  system. 

Let  x(n^,n2)  be  a  new  state  vector  defined  by 

x(nj^,n2)  =  T(n^,n2)x(n^,n2)  •  (5.21) 

Let  the  transformation  matrix  T(n^,n2)  be  nonsingular  for  all 
points  (nj^,n2)  in  the  region  S,  defined  in  (5.18)  ,  Making 
use  of  this  transformation,  equations  0.12)  to  (5.16). 
become 
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where 


x(n^,n2)  = 


x^(n^,n271 


X2(n^,n2) 


Xj^(n^,n2) 


(5.14) 


i  %• 


b(n^,n2)  = 


6^(n^,n2) 

S  2  *  ^2  ^ 


(ti^ ,  112) 


(5.15) 


^(n^^  1^2^  ”  1  * ^2 *^K "  (5.16) 


The  initial  conditions  for  the  model  above  are: 


x(n^,-l)  =  0  ,  x(-l,n2)  =  0  n^,n2  =  0,1,2, 


(5.17) 


If  the  input  array  {u(n2^  ,n2)}  has  a  finite  support 

S  =  { (nj^,n2)  ,  0  £  n^  <_  N^,  0  £  n2  £  N2}  ,  (5.18) 

and  we  are  interested  in  computing  the  output  y(nj|^,n2)  on  S, 
then  we  only  require  a  finite  set  of  boimdary  conditions , 
which  is  given  by 


x(nj^,-l)  =  0  for  n^^  =  0,1,.,., ,  (5.19) 


Figure  5 . 


Decomposition  of  rectangular  grid  for  obtaining 
recursion 


Figure  5,2:  Neighbors  of  (n, ,n2)  required  to  compute 
state  at  (n, .n^/  by  the  state-space  model 
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Let  us  consider  x^(nj^,n2)  for  an  arbitrary  but  fixed  i, 

i=l,2 ,  .  ,K.  Splitting  the  sumnation  in  (5.8)  according 

to  the  masks  in  Figure  5.1,  xve  have 

ni^-1  n2 

x.(nj^,n2)  =  ^  ^  3^(m^,m2)u(m^,m2)  + 

11I2 

ni  n2-l 

r  z  (5.9) 

m^=0  m2=0 

ni"!  n2-l 

-  Z  Z  6^(m^,m2)u(m^,m2)  + 

m^=0  m2=0 

Bi(ni,n2)u(ni,n2) • 

Substituting  (5.8)  in  (5.9)  and  making  use  of  (5.7)  ,  we 
have 

x^(nj^,n2)  =  x^(nj^-l,n2)  +  x^(n^,n2-l)  -  x^ (n^-1  ,n2-l) 

+  g^(n^,n2)u(nj^,n2)  i=l,2 . K  .  (5.10) 

K 

y(nj^,n2)  =  a^(n^,n2)x^(nj^ ,n2)  •  (5.11) 

Equations  (5.10)  and  (5.11)  can  be  rewritten  as 

x(n^,n2)  =  x(n^-l,n2)  +  x(n^,n2-l)  -  x(n^-l,n2-l) 

(5. 12) 

+b(n^,n2)u(n^,n2)  > 


y(n2^,n2)  c(n^,n2)x(n2^,n2) 


(5.13) 


Let  h(n2^,n2;  m^,in2)  impulse  response  of  a  2-D 

first  q-uadrant  quarter-plane  causal  discrete  LSV  system, 
and  let  u(n^,n2)  be  the  input  to  the  system  with  support  in 
the  first  quadrant.  Then, 

h(n^,n2;  m|^,m2)  =  0  for  n^  <  m^  or  n2  <  m2  ,  (5.3) 


u(n^,n2)  “  0  for  n^  <  0  or  n2  <  0  . 


(5.4) 


Substituting  (5.3)  and  (5.4)  in  (5.I)  the  input-output 
relation  is  given  by 

ni  n2 

y(n^,n2)  =  Z  Z  h(n^,n2;  m^,m2)u(m^,m2)  .  (5.5) 

m^®0  m2*0 


If  the  impulse  response  h(nj^,n2;  m^,m2)  is  a  K-th  order 
degenerate  sequence  given  by  (5.2),  then,  substituting 
(5.2)  in  (5,5)  we  have 

ni  n2  ^ 

y(n^,n2)  =  l  z  (  E  a^(n^,n2)  S^(m^  ,1112)  )u(m^  ,1112),  (5.6) 

m2^*0  m2“0  i®! 

After  some  manipulations  involving  the  interchanging  of 
summations ,  we  obtain 

K 

y(n^,n2)  =■  z^  a^(n^,n2)x^(n^,n2)  .  (5.7) 


where 


ni  n2 


x^(nj_,n2)  =  Z  Z  3^(m,  ,m2)u(m^,m2)  .  (5.8; 
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Linear  Shift- Invariant  systems  the  sum  in  (5.1)  gets  trans¬ 
formed  into  the  standard  2-D  convolution  and  there  are 
very  efficient  techniques  for  calculating  the  2-D  discrete 
convolution.  These  techniques  include  transform  domain 
methods  that  make  use  of  fast  algorithms  and  there  are 
also  state-space  models  for  recursive  implementation  in  the 
spatial  domain  [ch.  4  in  5.3]. 

For  LSV  systems;  the  transform  domain  techniques  cannot 
be  used  in  general  and,  in  addition  to  this,  there  are, 
as  yet  no  state-space  recursive  models  for  2-D  LSV  systems. 
Some  very  restrictive  1-D  state- space  models  for  LSV 
systems  have  appeared  in  the  literature.  They  permit  the 
analysis  of  1-D  LSV  systems  tmder  very  strong  restrictions 
imposed  on  the  impulse  response  [5.16'],  [5.1/]. 


5.2  2-D  Quarter  Plane  Causal  State-Space  Model 


In  order  to  develop  a  state-space  model,  the  following 
definition  has  to  be  introduced  as  a  natural  extension 
of  the  1-D  case  [5.14]. 

Definition  5.1:  A  sequence  {h(n^,n2;  m^,m2)}  is  a  K-th 
order  degenerate  sequence  if  it  can  be  expressed  as 


h(n^,n2;  ti^,m2) 


K 

Z^a^(nj_,n2)  3^(m^,m2)  , 


(5.2) 


where  a^(n^,n2)  (i»l,2 , . . . ,K)  and  0^(m^,m2)  (i=l ,2 , . . . ,K) 
are,  respectively,  linearly  independent  ftinctions  of 
(n^,n2)  and  (m^,m2) . 
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(where  the  impulse  response  sequence  was  assumed  to  be  either 
expressible  or  approximated  in  k-th  order  degenerate  form) 
domains.  In  thi«  chapter  attention  will  be  focused  on  the 
development  of  a  state- space  model  for  a  k-D  LSV  system  whose 
impulse  response  is  expressible  in  K-th  order  degenerate  form. 
For  convenience  in  exposition,  the  development  of  the  model 
will  be  initiated  for  the  2-D  quarter-plane  causal  case  and  a 
complete  proof  for  the  feasibility  of  extension  to  the  k-D 
(k  >  2)  case  based  on  a  (2^-1) -point  recursion  will  be  given. 
The  possibility  of  generalizing  the  model  to  cover  multi¬ 
dimensional  weakly  causal  LSV  systems  will  be  substantiated. 
The  model  for  the  inverse  system  will  be  derived  and  applied 
to  problems  in  deconvolution.  Nontrivial  physically  motivated 
examples  will  be  included. 

A  two  dimensional  discrete  LSV  system  can  be  described 
by  the  discrete  sum 

os  « 

y(n^.n2)  =  z  Z  h(n^,n2;  m^,m2)u(m^,m2),  (5.1) 

•  00  ’•00 

where  u(m)  is  the  input  at  coordinate  point  ra“(m^,m2), 
y(n)  is  the  output  at  coordinate  point  n=‘(n^,n2)  and 
h(n^,n2;  m^,m2)  is  the  response  of  the  sytem  at  the  point  n 
to  a  unit  impulse  at  the  point  m. 

The  input-output  relation  in  (5 . i)  is  far  from  being 
very  convenient  because  it  is  not  recursive  and,  therefore, 
the  amount  of  computational  time  required  for  implementing 
the  relation  will  be  very  large.  In  the  case  of  2-D  discrete 


sufficiently  general,  poses  some  dimensionality  problems  in 
implementation,  in  spite  of  the  fact  that  the  restorations 
besides  being  positive  are  spatially  smooth  and  not  overly 
sensitive  to  noise  in  the  image  data  at  object  points  that  are 
near  the  background  radiance.  Then,  there  is  the  class  of 
iterative  methods  based  on  gradient  type  [5.9]  and  conjugate 
gradient  type  [5.10]  of  optimization  procedures.  Iterative 
methods  are  not  very  suited  for  image  restoration  because  of 
the  problem  of  noise  suppression.  If  the  blurred  image  to 
be  processed  includes  a  noise,  then  this  noise  strongly 
deteriorates  the  quality  of  the  restored  image  as  the  number 
of  iterations  increase.  Thus,  a  method  is  desired  that  enables 
one  to  suppress  noise  amplification  more  efficiently,  while 
restoring  the  image  sharply.  The  three  ways  to  suppress  noise 
amplification  are:  (1)  stop  the  iterations  at  a  moderate 
iteration  number;  (2)  introduce  constraints;  and  (3)  reblur 
the  blurred  image  as  done  in  [5.11]  for  LSI  systems.  In  spite 
of  these  efforts,  the  general  problem  of  space- time  computa¬ 
tional  complexity  has  limited  the  use  of  practically  all 
methods  (including  those  already  cited  and  others  [5.12], 
[5.13]),  in  shift-variant  multidimensional  problems. 

A  more  general  approach  than  that  permitted  via  use  of 
coordinate  transformation  was  pursued  recently  in  [5.14], 
where  the  analysis  was  restricted  to  one- dimensional  (1-D) 
discrete  LSV  systems.  Investigations,  however,  were  carried 
out  both  in  the  frequency  (using  the  discrete  version  of 
Zadeh's  generalized  transfer  function  introduced  in  the 
discussion  of  continuous  time  systems  [5.15])  and  time 


The  last  row  of  B.  (n, )  in  (5.40)  can  be  expressed  as 

[3.(0)  3.(1)  ...  8.(n3_)]  =  [8.(0)  3.(1)  ...  3.(n3_-l)0]  + 

[0  0  ...  0  l]3^(n^)  .  (5.41) 

Making  use  of  Definitions  5.3  and  5.4.  the  expression  in  (5.41) 

can  be  written  as 

0  0 

B.  i(ni)  =  ^1(1)  B.  i(nT)  +  ^.(nT) 


(5.42) 
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which  corresponds  to  the  expression  in  (5.39)  for  k=l. 

From  the  matrices  in  (5.34)  to  (5.36) ,  for  any  integer  k  ^  1, 
we  can  write 


[0  0...g  (n2.n3....n^^p] 


0roi.3...i.,.i  lx.^x.>...  1.  , -1  ^ 

_  f,  2  3  k+1  ^  2  3  k+l.p  r, 

^^+1  “  hi+1  ^-i,k+l^’^l’^2’ •  •  • ’’^k+l)  • 


If  we  now  assume  that  the  hypothesis  is  valid  for  a  fixed 
integer  k,  with  n^  fixed,  we  can  write 


00.  .  .0 

^k 


2i:k<^2 


•  ’"k+1^  “ 


^  i+1 


'^k+1^ 


00. . .0 


(5.44) 


Also  from  (5,43),  after  setting  i2”^3” ' ‘ ’ ~^k+l*^ ’  have 

00. ..0  100. ..0 


^+1  -i.k+l^’^l . ^k+1^  ^+1 


-i,k+l^"l . ^k+1^ 


0  0  ...  0  n. 


+  [0  0  ...  g  ^ 


-i,k^^2 . ’^k+l^  ^  • 

(5.45) 


Substituting  (5.44)  in  (5.45),  we  have 


00. . .0 


10. . .0 


^+1  -i,k+l^^l' •  •  • ’^k+1^  ^+1  -i,k+l^^l . ’^k+l^ 


+  ^^(-Dj'^^ig  0..  .c^(j)B^^j^(n2 . n^^^^)  ; 


00.  .  .0 


+  ig  g. .  .g 


^k+l^’^l . ’^k+l^  ^  • 


(5.46) 
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From  (5.43),  making  use  of  Definitions  5,4  and  5,5,  we 
obtain 


to  0...0  2^(J)  (nj . 


0  1 

"  24c+l^j^^^i,k+l^^l’ •  ' ’’^k+l)  .  (5.47) 


Also,  from  Definition  5.3,  we  have 


00.  .0 

[0  0.  .  .0  ^k+l^’^l’  •  ‘  '^k+1^  •  (5.^8) 


Substituting  (5.47)  and  (5.48)  in  (5.46),  we  have 

00.  .0 

^+1  -i.k+l^’^l' •  ”^k+l^ 

10..  0  k  i+i  0  1 

^•^+1  ^-i,k+l^’^l' •  • ’’^k+l^ 

00.  .0 


■'■  -k+1  -i,k+l^’^l’  •  '  ’’^k+l^ 


(5.49) 


From  Definition  5.4,  it  is  straight  forward  to  show  that 


10. .0  0 

+  o_j^+l(l).  j=l 


-k+1 


0  1 

H4c+l<j)  =<  ^k+l^j^  ik+l^j-^^*  j=2,3 . k  (5.50) 


0^+1  (k),  j=k+l 


Consider  now  the  expression  in  square  brackets  in  (5.49) 
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which  corresponds  to  the  hypothesis  for  k+1.  Therefore, 
by  the  principle  of  mathematical  induction,  the  proof  is 
now  complete. 

The  following  definition  has  to  be  introduced  at  this 


point: 


Definition  5.6:  Define, 


(^) 


' ^2  ’  ■  ■  ■  ’ ^  '^i  ^  ’  ^2 "^s  ^ ^^  ’  ■  ■  ■  ’ ^k” ^s  ^ 

s*l 


where  z  r  (i)=j,  r  (i)=0  or  1  for  each  i  and  the  sets 
i=l  ®  ^ 

{r^(l),  r^(2),....  r„(k)},  s=l,  2 ,  .  .  .  ,  (^)  are  mutually  distinct. 

S  S  w  J 


Theorem  5.2:  Given  a  k-D,(k  ^  1),  positive  cone  causal  LSV 
system  whose  impulse  response  is  a  K-th  order  degenerate 
sequence  characterized  by  expression  (5.30).  Then,  the 
superposition  sum  in  (5.29)  can  be  implemented  by  a  K-th 
order  state-space  model.  This  state-space  model  is  based 
on  a  (2^-1) -point  recursion  described  by, 


x^(n^,..,n^)  =>  z^(-l)^'^^xi(n]_, .  .  ,n|^)  +  3^(n^,  . .  ,n^)u(n^,  .  .  ,n^) 

. K,  (5.54) 


y(n^, . . ,n^) 


K 

.i:  ai(ni,...nk)Xi(ni,..,n^). 


(5.55) 


Proof:  From  expressions  (5.33)  to  (5.36)  ,  making  use  of 
Definition  5.3,  we  have 


x^  (n^^- i 2^ ,  n^— i^  >  • 


■“k-V 


^1^2 ’ ’ ^k 

^  (5.56) 

i«l,2 . K  . 


Substituting  the  expression  (5.39)  of  Theorem  5  . 1  in  (5.56), 
after  setting  i^=i2“.  .  .  *i^=*0  , 
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Xi(ni,  .  .  ,Xj^)  =  ^ 

00.  .0 


6i(ni....nj^)Uk  . 


(5.57) 


From  (5.56),  making  use  of  Definitions  5.k  and.  5. 6,  we  have 


x^(n^,...n^)  *  -i,k^^l’ •  ' •  (5.58) 


From  the  ordering  defined  for  the  vector  we  have 

00.  .0 

“(“r-'-V  ■  ik  • 


(5.59) 


Substituting  (5.58)  arid  (5.59)  in  (5.57),  we  obtain  the 
expression  in  (5,54)  The  recursion  is  based  on  i-points, 
where 

k 


n  -  E  (^)  -  2^-1  . 
j-1  J 


(5.60) 


The  proof  is  now  complete. 
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The  example  below  serves  to  clarify  the  various  notations 
introduced  in  the  preceding  discussion.  For  brevity  in 
exposition,  this  example  tackles  the  k=2  case,  for  which  the 
3-point  recursion  arrived  at  in  section  5.2  is  verified  as 
a  special  case  of  the  general  result. 

Example  5.1  Consider  (5.32)  for  the  case,  k=2 .  Then,  the 
expression 

ni  n2 

x^(n^,n2)  =  E  E  B^(m^,m2)u(m^,m2) , 
m^=0  m2=0 


i 


is  considered  for  the  case  when  n^=n2=2.  Then  counterpart  of 
the  matrix  form  representation  in  (5.33)  is  given  below  in 
expanded  form.  Note  that  the  lexicographical  ordering  described 
has  been  adopted.  It  is  clear  that  in  the  k=2  case,  this 
ordering  results  from  a  row-by-row  scan  of  each  of  the  arrays 
{Xi(ii,i2)}  and  {u(i^,i2)}. 


I  X  (0,0) 

i  ' 

’  !  X.(0,1) 

i  ^ 

!  X  (0,2) 

I  ^ 

i  x(l,0) 

1  I 

I  x.(l,l) 

;  X  (1,2) 

I 

;  X  (2,0) 

I  ' 

'  x.(2,l) 

■  x.(2,2) 

I  i  I 


3^(0, 0) 

3^(0, 0)  8. (0,1) 
8^(0, 0)  3^(0, 1) 
S.(0,0) 

3^(0, 0)  3^(0, 1) 
'  3. (0,0)  3. (0,1) 

I  ^ 

!  3. (0,0) 

3^(0, 0)  3. (0,1) 
3. (0,0)  3. (0,1) 


8. (0,2) 

3^(1, 0) 
3^(1, 0) 
3^(0, 2)  3. (1,0) 
3. (1,0) 
3j^(l,0) 
3. (0,2)  3. (1,0) 


6. (1,1) 

3.(1, 1)  3. (1,2) 

3. (1,1) 

8. (1,1)  3. (1,2) 


3. (2,0) 

3. (2,0)  3. (2,1) 
3. (2,0)  3. (2,1) 


u(0,0) 
u(0,l) 
u(0,2) 
u(l,0) 
u(l,l) 
u(l,2) 
u(2,0) 
u(2,l) 
3. (2,2)  u(2,2) 


(5.61a) 


I 
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Specializing  the  representation  in  (5.33)  to  this  case, 
(5.61a)  can  be  written  as, 

X2  =  2’ 


where , 


and 


Bi  i(n2) 

-i,1^^2^ 

i(n2) 

B°  i(n2) 

b1  i(n2) 

Si(j.O) 

6i(j,0) 

3^(3.!) 

^^(j.O) 

Sj_(j,l)  . 

Bi,i(n2) 


6^(j ,n2) 


(5.61b) 


(5.61c) 


Clearly,  (5.61b)  and  (5.61c)  give  the  relevant  specializations 
of  the  matrices  in  (5 . 34) - (5 . 36) .  Applying  Definitions  5.3 
and  5.4  to  the  case  under  consideration,  we  have 

00 


±2  =  [000 

000 

001], 

£2°  =  [000 

001 

000] , 

5°^  =  [000 

000 

010], 

4^  =  [000 

010 

000], 

a_2(l)  =  ^2^  ^2^ 

£2(2)  =  . 
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Let  r.  denote  the  j-th  row  in  the  matrix  B.  ^  of  (5.61a), 
j=l,2, ... ,9.  Clearly 

I9  =  ES  ^  ■  -5  °  3^(2,2)  ] 

the  above  equation  can  be  rewritten  in  the  form, 

-2°  2 ^’^1*^2^  ^  -2^  -i,2^’^l’’^2^  , 2 ^’^l ’^2^ 


-  5^^  i2°  2.  (2,2) 


+  ^2°  3^(2, 2)  , 


(5.62) 


which  is  the  relevant  specialization  of  (5.39). 

Note  that, 

Xi(ni-i,n2-j)  =  £2^  B^  2  >’^2^  ‘ -2  ’  (5.63) 

Multiplying  (5.62)  from  the  right  by  U2  and  using  (5.63), 
we  get 

x^(n^,n2)  =  x^(nj^,n2-l)  +  x^(n^-l,n2) 

(5.64) 

-  x^(n^-l,n2-l)  +  6^(nj^,n2)u(n^,n2). 
Applying  Definition  5.6, 

xl(n^,n2)  =  x^(n^-l,n2)  +  x^(n^,n2-l),  (5.65a) 

x?(n^,n2)  x^(nj^-l  ,n2-l)  .  (5.65b) 

Substituting  (5.65a)  and  (5.65b)  in  (5.64)  we  get 

1  2 

x^(n^,n2)  =  x^(n^,n2)  -  x^(nj^,n2)  +  ^  (n^  ,  02)  u(n^  ,  02)  , 

which  is  the  relevant  specialization  of  (5.54). 
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5.3.2  Extension  of  the  k-D  State-Space  Model  to  a  Causality 
Hypercdne 

The  state-space  model  developed  in  Section  5.3.1  for 
a  k-dimensional  positive  cone  (or  hypercone)  causal  discrete 
LSV  system,  can  be  naturally  extended  in  order  to  represent 
a  k-dimensional  (k  >  2)  weakly  causal  discrete  LSV  system. 
Definitions  5.7  and  5.8  given  next,  are  introduced  in  order 
to  reach  the  desired  goal. 

Definition  5.7:  A  k-dimensional  causality  hypercone  C^;  is 
the  intersection  of  k  half  hyperplanes  Hp^^,p^2' • • ’Pik 
(i»l, 2 , . . . ,k) ,  where 


'^Pil' Pi2' ■  ■  '  ' Pik  *  ( .  .  .  ,Xj^)  I  (x^,  .  .  .  ,Xj^)  e  R  , 


Pi  A  Pik=ht  i 

(1-1,2 . k) 


and  p^j(i,j»l,2 . k)  are  non-negative  integers, 

satisfying 

det  «  »  1, 


where 


Pll  Pi2  • • •  Plk  ! 

f 

A  ^ 

y  ^ 

Pkl  Pk2  •  •  •  Pkk 


(5.66) 


It  is  important  to  note  that  in  a  k-dimensional  causal¬ 
ity  cone  C^,  any  vector  v  going  from  the  origin  to  a  point 
P  in  can  be  expressed  as  a  linear  combination  of  the 
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vectors  i=l,2,...,k.  These  vectors  are  called  the 

generator  of  the  causality  cone  and  it  can  be  shown  that 
the  generator  , i=l, 2 , . . . ,k  corresponds  to  the  i-th 
column  of  the  adjoint  of  the  unimodular  transformation 
matrix  i . 

Definition  5.8:  A  k-D  discrete  LSV  system  with  impulse  , 

response  h(n^, . . . ,n^;m^ . causal  on  a  causality 

hypercone  if  and  only  if 

h , n^ ;  m^ ,  -  0  for  ( 2, ^ ,  .  .  .  ,  i  ’ 

with  2^  =  n^  -  m^  (i=l,2 . k). 

The  one-to-one  and  onto  mapping  « ,  defined  in  (5.66), 
maps  any  integer  point  in  C  onto  a  unique  integer  point  in 
the  positive  cone 

n  ^  Qt  n  z^  and  $  [  (0, 0  .  .  .  ,  0)  ]  =  (0 . 0 .  .  .  ,0) , 

c  i 

Given  a  k-D  discrete  LSV  system  with  a  K-th  order  degenerate 
impulse  response,  which  is  causal  in  C  ,  then  use  of  the 

mapping  $  defined  in  (5.66)  inverse,  a  state-space 

model  can  be  derived.  The  resulting  K-th  order  state-space 

model  will  be  based  on  a  (2  -1) -points  recursion.  For  the 
sake  of  clarity  and  brevity  in  exposition,  the  procedure  is 
described  for  the  2-D  case,  from  which  the  k-D  (k  >  2) 
counterpart  can  be  obtained  as  a  direct  generalization. 
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y(nj^,n2)  “  l  Z  h(n,,n2;  m^,m2)u(m^,m2) 

(ni^^  I  ni2 )  £ 

(ill, 1X2)  £C^ 


C5.67) 


be  the  input-output  relation  of  a  discrete  LSV  system,  and 
let  ^h(nj^,n2 ;  mj^,m2)}  be  causal  in  C^.  In  addition  to  this, 
assume  that  ih(n^,n2 :  m^,m2)}  is  a  K-th  order  degenerate 
seqxience  expressed  as  in  (5.2).  Using  the  map  $  in  (5.66) 
with  k=»2,  we  map  the  input,  the  output  and  the 
impulse  response  in  (5.67),  as  follows 


u(m^,m2)  *  u(m^,m2) 


“1 


®1 


(5.68a) 


y(n^,n2)  “  y(n^,n2) 


^1 

^2 


<4) 


(5.68b) 


ia(n^,n2;m^,m2)  = 


h(n^,n2 ;m^,m2) 


®1 

’^1 


-1 

-1  ^1 
^2 


(5 . 68c) 


Then,  h(n^,n2;  mj^,m2)  is  first  quadrant  quarter  plane 
causal.  The  input  array  u(m^,m2)  and  the  output  array 
y(n2^,n2)  have  support  on  Q^.  In  addition  to  this,  from 


5.30 


(5.2)  and  (5.68c)  h(nj^,n2  i-s  a  K-th  order  degenerate 

sequence  ,  and  it  can  be  written  as 


K 

a^(n^,n2)  3^(m^,ni2)  , 


(5.69) 


where 

ct^(n^,n2)  ~  ^^(n2^»n2)  j 

V"*  V 

3^(m^,m2)  =  3^(m^,iii2) 


From  (5.67)  and  (5.68),  we  have 

^2  .  . 

y(n2_,n2)  =*1  Z  h(nj^,n2;  m^,m2)u(m^ ,1112)  .  (5.71) 

mi“0  m2“0 

For  the  first  quadrant  qioarter  plane  causal  LSV  system  in 
(5.71).  with  its  impulse  response  given  by  (5.69),  we  can 
now  write,  using  equations  (5.12)  Co  (5.15).  ^  state- 
space  model  with  support  on  Q^. 

x(nj^,n2)  *  x(n^-l,n2)  +  x(n^,n2-l) 

(5.72) 

-  x(n^-l,n2-l)  +  b(n^,n2)u(n^,n2)  . 
y(n^,n2)  =  c(n^,h2)x(n^  ,112)  . 
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where 


b(nj^,n2)  = 


3  ^  ^ 


3  >  ^2  ^ 


(5.73a) 


’  ^2  ^  ^K  C'  7  3b ) 


The  initial  conditions  are 


x(n^,-l)  =  0,  x(-l,n2)  =  0  n^,n2  =  0,1,2,... 


(5.73c) 


Mapping  back  to  the  state-space  equation  described  in  (5.72) 
and  (5.73),  (for  notational  brevity,  replace  P]_i  >  P]_2  *  ^21  ’  ^22 
in  (5.66)  by  p,  r,  q,  t,  respectively) 

x(n^,n2)  =  x(n^-t,n2+q)  +  x(n^-fr  ,n2-p) 

(5.74a) 

-  x(n^+r-t,n2+q-p)  +  b (n^ ,n2)u(nj_ ,n2)  ,  (n^,n2)  e  . 

y(n^,n2)  =  c(n^,n2)x(n^,n2) ,  (n^,n2)  e  ,  (5.74b) 


where 


b(n^,n2)  = 


S^(n^,n2) 


(5.75a) 


» *■‘2' 


•  •  ?  ^^(u2^,n2)  ] 


(5.75b) 


The  initial  conditions  are 


x(tn+r , -qn-p)  =*  0,  x(-m-t,pn+q)  =  0 
n®0 , 1 , . . . 


(5.75c) 
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=  -f-  ‘'o'— 


r . cos (9 . -9^) -r 


r^sin(0,-9„) 


-),  (5.90) 


where  the  function  hQ(x^,y^)  is  the  response  to  an  impulse 
at  Thfi  form  of  hgC*)  for  several  special  cases 

such  as  spherical  and  coma  aberrations  are  derived  in  [5.20] 
For  Cartesian  coordinates  (x^,X2)  and  of  the  image 

and  object  planes,  respectively,  we  have 


x^  =  r^coso^ 

T,  =  r^cose 
1  0 

^2  “  r^sine^ 

To  =  r^sine 
L  0 

(5.91) 


Substituting  (5.91)  in  (5.90),  the  impulse 

response  in  Cartesian  coordinates,  h(x^ ,X2 : , T2) ,  for 


the  coma  aberration  becomes 


■1  XnT.,^XrtTrt  XoTt  **X'nT^ 

h(x^,X2;Tj_,T2)  -  -2— -■-j  hQ(^4^ - (5.92) 

Ti  +  t2  +  T2  +  t2 


where  hQ(x^,X2)  is  the  response  to  an  impulse  at 
that  is,  hQ(X|^,X2)  *  h(x2^,X2;1.0)  ,  is  given  by 


hQ(x^,X2)  = 


'x^  -  3X2 


'x^  -  3x2 


for  (x^,X2)  e  I  , 


for  (x^,X2)  £  II  , 


(5.93) 


derived  in  (5 . 88)  -  (5 . 89)  ,  a  3tat-?-bpace  model  of  the 
inverse  system  was  implementec. 

The  (32x32) -points  original  object,  shown  in  Figure 
(5.4a)  was  blurred  by  the  simulated  motion  blur  described 
in  (5.86)  with  the  parameters  in  (5.82)  and  (5.84)  taken 
typically  as  T=l,  a=2,  a=l  and  i=.2.  The  resulting  motion 
blurred  image  is  shown  in  Figure  (5.4b).  '  From  the  motion 
blurred  image  in  Figure  (5. 4b)  by  using  a  (2-K+l)-th,  (K=10)  , 
state-space  model  of  the  inverse  system,  the  original  object 
was  reconstructed.  The  reconstructed  object  is  shown  in 
Figure  (5.4c),  From  Figure  (5.4c),  we  can  conclude  that 
for  a  fairly  low  order  state-space  model,  the  reconstruction 
is  very  accurate , 

5.5.2  Coma  Aberration 

Within  the  geometrical  optics  model  of  image  formation, 

aberrations  are  described  by  the  ray  aberration  function 

which  is  a  vector  in  the  image  (output)  plane  from  the 

Gaussian  image  point  to  the  point  where  the  ray  actually 

intersects  the  image  plane.  The  Gaussian  image  point  is 

the  image  in  the  output  plane  of  a  point  in  the  input  plane 
in  absence  of  optical  aberrations.  Let  (r.,e.)  and  (r  ,e  ) 

denote,  respectively,  the  polar  coordinates  of  the  image 

and  object  planes.  It  has  been  shown  in  a  paper  by  Robbins 

and  Huang  [5.19]  that  the  impulse  response,  h(r^ , 9^ : r^ , e^) , 

for  coma  aberration  is  given  by 
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RBaaaaaaeBaaaaaBaaaeaaeaaaBaeeaaeBaaBBaaaaaBeaaaaaaaaBaaaaaaaaaa 

BBaaaaaaaaBeaaaeaaBaaaaaaaRaaBBeaeaaaBaeaaaaaeaaaaaaaaeaaeaaaaaa 

BBaaaaBaBaaaeaaaaaaaeBaEeaeaaaaasaeeasaaaaaaaaaaaaeBaaaeaeaaaaaa 

RaaaaaeaaeBeaaBBeaaBBRaBaeaBeaaaaeaaBaaeBBaaBeaeaBBaBaaaaaBaaaaa 

RBaaaaaaaa  - - — . . --aesaRBEEaaaeHHBaeeae— . . : :  — . . — BBaaaaaaaa 
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Figure  5.4c:  Restored  object  from  motion  blurred  image 
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iitefSEZfeiefBicesfCffeeiiBSBXffeifSBffsxfSsssssffsxssssssssssss 
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BBSSSf  ff8BBBBffff??ffffHH99^f?^CC*''*’fr.t«f-f  ;  ;  ,  .  _ ,  .  !  --f  + 

BB888S8fiB8S8SSBBHH99Mf^9C'**^t1^*h'f  t  t  .  «  _ ,  •  t  t  4 

B8BB8B988SSS18tf<MtHH99MM0C*^t:^4*f  t  t  «  •  t  !~~44 
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BIBBBBBXXBBXXBIBXBEEBBBBBXeBBBBXBBBBXXEBXXiBeXBBBXBBBXXZBBXXXXXX 
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Figure  5.4b;  Motion  blurred  image 
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Figure  5.4a;  Original  Object 
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(NxN)  of  the  input  array.  The  discrete  motion  blurring  is 
carried  out,  line  by  line,  by  means  of  the  following 
expression 

^1  , 

y(n^,n2)  =  l  hjj(n]^,m]^)u(m^,n2)  .  (5.86) 

m^=0 

n^=0,l . N-1,  i=l,2 

Analogously,  before  implementing  a  discrete  state-space 
model  for  the  inverse  system,  the  approximated  impulse 
response  h*(x,u)  in  (5.84),  has  to  be  discretized,  and  the 
resulting  discrete  impulse  response  hg(n,m)  is  given  by 

hg(n,m)  =  h*(n-A,m-A)  ,  (5.87) 

where  A  has  been  previously  defined.  Then,  from  (5.87)  and 
(5 -84) - (5 . 35) ,  we  have 

2-K 

h*(n,m)  =  Z  a.(n)B.(m),  (5.88) 

^  i=0  ^  ^ 

where 

cxQ(n)  =  1  ,  £Q(m)  =*  aQ(mA)  ,  (5.89a) 

oi  .(n)  =  cos(n-^  i-A)  .  6 .  (m)  =  a.  (mA)  ,  i=l,  .  .  .  ,K  . 

Tp  (5.89b) 

a^(n)  =  sin(n-  y~  >  e^(m)  =  b^(mA)  ,  i*K+l . 2K  , 

P  (5.89c) 

From  the  (2-K+l)-th  order  degenerate  sequence,  {h^(n,m) } 
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b^(u) 


where 


(u+a) 


{a^cos(w^(u2-d2))cos(w^d)  +  cos(w^u^)} 


_1_  ,  ^ 
D2  ^  T 


(5.35c) 


a 

(u+a) 


{S^cos(w^(u2-d))cos(w^d)  +  ^  cos  (w^u^) 


d  .  ,  TT  .  1  _  2k 

“  s  in  \  /  t  ^ 

k=1.2,...  . 


Wk  -p  . 


2-T(Wk  +  D^) 


1  N 


D2/2  . 


In  order  to  simulate  the  motion  blurring  phenomenon  on 
a  digital  computer,  the  impulse  response  h(x,u)  in  (5.83) 
was  discretized  and  hQ(n,m)  corresponds  to  the  discrete 
version  of  (5.83),  where 


hj^(n,m)  =  h(nA,  mA)  , 


for  a  fixed  A,  such  that  U  =  (N-1)"A.  The  constant  U 
has  been  previously  defined  and  N  corresponds  to  the  size 
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where 


Tp  >  Uj  -  U  +  D^, 


=  I±a  u  +  ^ 


and  will  be  defined  to  be  Tp/4.  Under  this  assumption 

h(x,u)  can  be  expanded  in  a  Fourier  series,  which  is 
truncated  to  yield, 


•  _  /2it 


h*(x,u)  =  ar,(u)  +  E  {a.  (u)cos(^kx) +b,  (u)  sinC^-  kx) } . 

k=l  -^p  p 

(5.84) 


The  a^(u)'s  and  b^(u)'s  in  (5-13),  obtained  in  closed  form, 


agCu)  = 


(U2  ■  •  “*2^' 


(5.85a) 


a^(u)  = 


a 

(u+oT 


(u+cx"5" 


{a^sin(w^(u2-d))cos(w^d)  -  sin(w^u^)  } , 


l_  a  ^ 

D2  ^  T  ' 


(5.85b) 


1 


{ 6j^sin(w^(u2-ci) )  cos  (w^d)  -  ^  sin(wj^u^) 


1  2k 


a  /  u  \  J-  _ 

■  T  ^2^  ’  D,  T 


k=l,2,...  , 


ulse  response  of  LSV  system  modeling 
blur 


over  these  superfluous  points.  The  discontinuity  at 

OL  T 

X  =  — —  u  +  ^  ,  will  be  slightly  smoothed  for  better 

Cl  a 

results  in  the  approximation.  From  the  above,  the  function 
to  be  approximated,  h(x,u) ,  is 


1 

a 

/  U+a 

u^  _<  X  £  ^2*^2  ’ 

h(x,u)  =  i 

1  sfc-  • 

^2“^2  —  ^  —  ^2  ’ 

(5.83 

1 

^  0 

elsewhere  , 

where 


u,  =  u  -  D, 


D2  =  (u2-u)/10 


w(x)  is  a  Hanning  window,  given  as 

w(x)  =  .5(1  -  cos(^  (X-U2)). 

The  constant  will  be  defined  later  in  this 
section.  The  function,  h(x,u),  is  plotted,  for  a  fixed 
but  arbitrary  value  of  u,  in  Figure  5.3. 

The  u  variable  will  take  values  over  a  finite  interval. 


[0,U];  then,  for  any  arbitrary  but  fixed  u  e  [0,U],  we  can 
consider  h(x,u)  to  be  periodic  in  x  for  a  sufficiently 
large  period,  T  . 
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5 . 5  Examples  of  Applications 

In  this  section,  two  examples  are  presented  which  illustrate 
the  application  of  the  state-space  model  for  the  inverse 
system  to  the  restoration  of  a  degraded  image,  where  the 
degradation  is  modeled  by  a  LSV  system  describable  by  (5.10) 
to  (5.17).  In  section  5.5.1,  the  physical  phenomenon 
responsible  for  image  degradation  is  motion  blur  while  in 
section  5.5.2  the  degradation  is  due  to  a  type  of  optical 
aberration,  referred  to  as  coma. 

5.5.1.  Space- Variant  Motion  Blur 

Consider  the  one- dimensional  motion  [5.1o] 


u  -  g(x;t)  »  a,a,u  >  0 

t  e  [O.T]  , 


(5.81) 


Using  Sawchuk's  analysis,  this  motion  is  modeled  by  a  LSV 
system  whose  impulse  response  is: 


h(x,u) 


d  T+a  1  oT 

u+a  ’  —  —  a  a 


0  elsewhere  , 


(5.S2) 


In  order  to  obtain  a  1-D  state-space 
model  for  this  blurring,  the  impulse  response  in  (5.82) 

vill  be  approximated  in  degenerate  form.  Before  carrying 

out  the  approximation,  it  is  important  to  note  that  the 

state-space  model  is  causal  and,  therefore,  it  will  not 

evaluate  the  impulse  response,  h(x,u) ,  at  points  x  <  u; 

this  will  be  used  to  our  advantage  in  extending  h(x,u) 


h(n^,n2;n^,n2).  0  for  all  0  ^  N  ,  i»l,2.  (5.79) 


5.4.2  Weakly  Causal  Case 

The  inverse  of  the  system  described  in (5 . 74)- (5 . 75) 
is  easily  shown  to  be: 

z(n^,n2)  =  A^(n^,n2)  [2(n^-t,n2-l-q)  +  z(n^+r,n2-p) 

-  £(n^+r-t,n2+q-p) ]  +  b^(n^,n2)y(n^,n2) ,  (5.30a) 

u(n^,n2)  =  c^(n^,n2) [z(n^-t,n2+q)  +  z(n^+r,n2-p) 

-  z(n^+r- t,n2+q-p) ]  +  d^(n^,n2)y(n^,n2) ,  .(5. 80b) 

with  initial  conditions 

z(tn+r , -qn-p)  =  0  ,  2(-m-t,pn+q)  =»  0  ,  n=»0,l,...  (5.80c) 

where  Aj(nj^,n2)  ,  bj(n^,n2)  £]^(nj_,n2)  and  d^(n^,n2)  have 
been  defined  in  (5.77). 

The  necessary  and  stifficient  condition  for  the 
existence  of  the  state-space  model,  in  (5.80)  is, 

Cj(n^,n2)  '  bj(n^,n2)  f  0  for  all  (n^,n2)  e  , 

or,  equivalently. 


h(n^,n2 ;n^,n2)  f  0  for  all  (n^,n2)  e  , 


z(n^,-l)  *0  ,  z(-l,n2)  =0  for  n^“0,l, . . . ,N  ,  i=l,2  ,  (5.76c) 
where 

Aj(n^,n2)  =  Ij^-b (n^ ,n2)  [c(n^,n2)b (n^ ,n2)  ] '^  c(n^,n2).  (5.77a) 
and  is  the  K-th  order  identity  matrix, 

bj(nj^,n2)  =  b(nj^,n2)  [c(nj^,n2)b(nj^,n2)  ]  (5.77b) 

c^(n^,n2)  =  -  (c(n^,n2)b(n^,n2) ]~^  c(n^,n2)  ,  (5.77c) 

d^(n^,n2)  =  [c(n^,n2)b(n^,n2) •  (5.77d) 

It  is  important  to  note  at  this  stage,  that  the  state- space 
model  of  the  inverse  system,  given  in  equations  (5.75)  And 
(5.77),  is  based  also  on  a  three  point  recursion.  Therefore, 
it  will  provide  a  very  efficient  deconvolution  procedure. 

In  addition  to  this,  for  a  reasonably  low  order  state-space 
model,  there  will  be  no  storage  problem  because  the  state- 
vector  is  only  two  dimensional  and  there  is  no  need  to 
store  it  over  the  entire  input  mask. 

The  state-space  model,  considered  here,  for  the 
inverse  system  exists  if  and  only  if 

c(n^,n2)b(nj^,n2)  0  for  0  ^  n^  1  N  ,  i=l,2,  (5.78) 

or  equivalently,  by  substituting  (5 , 15)  and  (,5 . 16)  in  (5.78) 
and  making  use  of  (5. 2)  ,  we  have  that  the  condition  in  (5.78) 
is  equivalent  to 
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From  the  above  equations,  we  can  see  that  the  state-space 
model  in  (5.74)  and(5_.75)  preserves  the  three  point  recur¬ 
sion.  It  is  important  to  emphasize  that  in  this  case 
the  recursion  in  (5.74a)  does  not  depend  on  the  closest 
past  neighbors  of  the  point  currently  under  consideration. 

5.4  State-Space  Model  for  the  Inverse  System 


In  section  5 . 4. 1,  the  state-space  model  for  the  inverse  of 
the  system  described  in  (5 . 10)- (5 . 17)  is  first  given.  In 
section  5.4.2,  the  coimterpart  in  the  weakly  causal  case  is 
considered. 


5.4.1  First  Qtiadrant  Quarter-Plane  Case 


Consider  the  state-space  model  of  the  first  quadrant 
qxoar ter -plane  discrete  LSV  system  described  in  (5 . 10)  -  (5 . 17)  . 
The  state-space  model  of  the  inverse  system,  whose  state- 
vector,  input  and  output  at  the  point  (nj^,n2)  are,  respectively, 
z(n^,n2),  y(n^,n2)  and  u(n^,n2)  is: 

^(n^,n2)  *  Aj(n^,n2) 

z(n^-l ,n2-l)  ]  +  bj(nj^,n2)y(n^,n2)  ,  (5.76a) 

u(n^,n2)  “  ^j(n^,n2)  [^(n^— l,n2)  ^(n2^,n2”l)  — 

z(n2^-l,n2"l)  1  +  dj(nj_,n2)y(nj^,n2)  .  (5.76b) 


with  initial  conditions 


!• 


where  C  is  a  constant  and  the  regions  I  and  II  are  clearly 
defined  in  Figure  5; 5. 

The  impulse  response  h(x^,X2;T ^,12)  for  the  coma 
aberration  in  (5.92)  is,  in  general,  non- causal,  but  by 
imposing  some  constraint  on  the  radius  Rq  in  the  pattern 
of  the  impulse  response  hQ(x^,X2)  in  (5.93)  as  shown  in 
Figure  5.5,  it  is  possible  to  perform  the  deconvolution, 
recursively.  The  restriction  on  Rq  corresponds  to  a 
limitation  in  the  amount  of  coma  aberration  to  be  tolerated, 
which  can  be  made  very  small  depending  on  the  quality  of 
the  lens . 

In  order  to  perform  the  deconvolution,  the  input  plane 
will  be  divided  into  its  four  quadrants  and  the  deconvolution 
will  be  carried  out  for  each  one  of  them.  The  final 
reconstructed  object  is  obtained  by  superimposing  the  four 
resulting  arrays.  For  this  procedure  to  be  valid,  without 
loss  of  generality,  an  object  point  outside  should  not 
affect  an  image  point  on  Qj^.  This  condition  is  illustrated 
in  Figure  5 . 6,  and  it  will  be  used  to  determine  a  bovind  for 

^0  ■ 

From  Figure  5.6,  the  inequality  d.^  <  ly|  is  sufficient 
to  guarantee  that  the  object  point  (1^,12)  will  not  affect 
any  image  point  on  Q^.  Simple  calculations  enable  one  to 
gat  the  maximum  allowable  Rg  required  to  recursively  restore 
an  (NxN)  input  array  to  be 


Rr 


(j)^  (N-1)  -  2 


(5.94) 
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As  an  example,  a  (31x31) -point  original  object  u(n^,n2) , 

—  ’^i  —  which  is  the  counterpart  of  the 

(32x32) -point  original  object  in  Figure  5.4a  was  blurred  by 

the  coma  aberration.  For  this  object  size,  N=31,  and  Rp, 

-  '^max 

.065.  The  blurred  image  y(n^,n2) ,  -15  i  1  15 ,  i=l,2 

was  obtained  by  performing  the  summation, 

15  15 

y(n^,n2)  =  z  Z  h  (n^,n2;m^,m2)u(m^,m2) 

m^*-15  m2— “15 

-15  1  1  15,  i=l,2 

where 


hQ(n^,n2  :m^,m2)  =  h(n^A  ,n2A  •,m^A  ,m2A) 

-15  ^  n^^  ^  15,  -15  ^  m^  £  15, 
i=l,2 


and 


A  -  1/15 

The  resulting  blurred  image  is  shown  in  Figure 5^ 7.  In 
order  to  implement  the  state-space  model,  the  impulse  res¬ 
ponse,  hjj(n^,n2;m^,m2)  ,  was  exactly  represented  as  a 
degenerate  sequence  by  means  of  the  DFT  of  hjj(n^,n2  ;m^,m2) 
for  each  pair  (m^,m2),  -15  ^  m^^  ^  15.  This  procedure  is 
explained  in  [5.21].  The  symmetry  of  the  impulse 

response  of  the  coma  blurring  was  used  in  reducing  the 
amount  of  DFT's  required  for  this  representation. 

A  state-space  model  of  the  inverse  system  was  obtained  and 
used  in  the  deconvolution.  The  deconvolved  image  matched 
exactly  the  original  image. 


BBBBBBNNBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBNNBBBBBBBBBB 
BBaB  =  >IIIIHIWHtNIBBIBBBBBBBBBBaBBaBBBBaBI8flMWHIINX»I==BBBBBBBB 
aaBBXXHHKXXXWNBMBBBBBBaBBaBBBBaBBBBBaaMMWMKXXXHHXZBBBBBBBa 
BBBBMMxxMXHHBaaaaaaaaaaBBBBaaaBaaaBBaBaaaBHHxxxxMMBBBaaaBa 
BaBBMMXXMMaBHHaBBBaBBBBaaBBBaaaBBBBBBBBaHHBBMMZXMMBBaBBBBB 
BBBBXX)(KltlNiaaHIWMMaBBBBaBaBBBaBaBBBaBBaBH<HIWNIBBMHIKXXZBBBBBaBB 
aaBB: :8eHHBBBaBBBBaflBBBBBaMBBaBBaBBBBBBBBBBHHBe: tBBBBBBBB 
BaBBHM  XXMMIHHiaBaBBBBBaaBBBaBBBBaBBBBBBaiHIHMMKX  HHBBBBBBBB 
aaaazzhMOo-f-f — xxBsaaaaaaaaaaaaaaaaaaBBxx — ffooMMZzaaaBaaaa 
aBBBeeuviHMSsiieeflaaBBaBflaBaaaBBiBBBBaaeeiissMMuueeaaBBBBBB 
BaaBZZHMNNZZXKHHaaaBiaBaBBBBBaBBBBaaaaHHXKZZNNHHZZaBBaBBBB 
BBBBMMXXXXZZKXHHaBaBaaBBBflflaBBaBBBaBBBHHXXZZXXXMlHMBBaBBBBB 
BBBBZZHMNNZZXXHHflBBBBBBBBBaBBBBBBBflBBBHHKXZZNNnMZZBBBBBBBB 
BBBBaSUUMMSSIieeaBaBBBBBBaBBBBBBBBBBBBeeilSSMMUUeSBBBBBBBB 
8B8BZZHMQ0-t-t— XZBBBBZXZZSSIISSZZZXBBBBXX— f+OOKHZZBBBBBBaa 
BBBBMM  KXWMMM8BMRXXZZWII-(-fX»-^fWMZZZXMM8BN<INMMKX  hMBBBBBBBB 

bbbb:  tesHHBRaaiHHiBBMWHHxxuuuuwuxxHHMWBBiiwBBBRHHBe:  :bbbrbbbb 
■BBBXXXXHIMBRNIHieBMIHBBXZBeHHHHHHBBXXBBIHIHeBHIHIBBHIMXHXXBBBBBBBB 
aaBaMMXXMI«BaHHWMBBI«W--lIRXXM»(XII—BWBBMMHHBBWWZXWMBBBBBBBB 
RBBBMMXXXXHHBRBBBBMIll-l-i-SSZZZZZZSS-l-t’MMBBaRBBHHXXXXMMBBBBBBRB 
aaBaxXXXHHHHaaWMXXXXUUeBZZMWZZBBUUXXXXIIIWBBHHHHZZXXBBaBBBBa 
BBBBBBBBBBBBBBRBBBBBBBRBBRBBBBBBBBBBBBBRBBBBBBBBBBBBBBBBBB 
BRBBflBZZBBBIBBBBBBBBBBBBBRBBBBBBBBBBBBBRBBBBBBZZBBBBBRBIBB 
aaBBBBBBBBBBBBBBBBBBBBRBBBBBBaaBBBBBRBBBBBBBBBBBBBBBBBBBBB 


Figure  5.7;  Image  blurred  by  coma  aberration 


5.6  Conclusions 


Procedures  for  compensating  for  effects  which  degrade  the 
accuracy  of  remotely  sensed  datas  by  mathematically  inverting 
some  of  the  degrading  phenomenon  are  required  in  biomedical, 
industrial,  surveillance  and  earth  and  space  applications. 
Images  to  be  restored  are  often  degraded  by  a  linear  spatially 
varying  operation.  Degradations  due  to  motion  blurring  and 
optical  system  distortions  often  require  the  imaging  systems 
to  be  analyzed  by  modeling  the  degradation  as  a  linear  shift- 
variant  operation.  High  speed  digital  computers  have  been 
primarily  responsible  for  the  development  of  techniques 
for  image  restoration  of  shift- invariant  motion  blur,  but 
counterparts  of  such  techniques  in  the  shift-variant  case  are 
severely  handicapped  by  the  required  space- time  computational 
complexity.  Attention  has  been  directed  to  the  alleviation 
of  this  shortcoming  and  the  results  arrived  at  are,  first, 
briefly  summarized  and  then  directions  for  research  in  the 
immediate  future  are  provided. 

For  any  2-D  discrete  first  quadrant  quarter-plane  causal 
linear  shift-variant  (LSV)  system,  who$  impulse  response  is 
a  K-th  order  degenerate  sequence  a  K-th  order  state- space 
model  was  obtained.  This  model  is  recursive  and  is  based  on 
a  three- term  recurrence  formula  relating  any  point  in  the 
state-space  model  to  its  three  closest  past  points  and  the 
current  input.  The  state-space  model  was  extended  in  order 
to  model  2-D  discrete  LSV  systems  with  support  on  a  causality 
cone.  Subsequently,  the  2-D  quarter-plane  causal  and  weakly 
causal  discrete  models  were  extended  to  the  n-D  (n  >  2)  case. 
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The  resulting  state-space  models  are  recursive  and  based  on 
a  (2^-1) -points  recurrence  formula,  which  for  the  causal 
case  used  the  (2^-1) -closest  neighboring  past  points  in 
addition  to  the  input  in  order  to  compute  any  current  output 
state.  For  the  weakly  causal  case,  the  (2^-1)  computed 
points  required  are  not,  in  general,  the  closest  neighbors  to 
the  present  output,  which  is  being  computed.  Conditions  for 
the  existence  of  a  2-D  state-space  model  for  the  inverse 
system  are  readily  derivable  from  the  original  one.  Models 
for  the  2-D  LSV  system  and  its  inverse  can  be  used  to  perform 
analysis  and  deconvolution  problems  very  efficiently.  This 
can  be  substantiated  from  dervided  expressions  for  space- time 
computational  complexities  [5.21], 

Examples  of  physically  motivated  applications  making 
use  of  the  theoretical  results  developed  have  been  worked 
out.  These  applications  include  effects  of  1-D  LSV  motion 
blur  and  the  blurring  due  to  Seidel  aberrations  of  a  lens; 
in  particular,  the  2-D  LSV  coma  aberration  was  studied  in 
detail.  The  reconstruction  of  the  original  object  from 
the  LSV  blurred  image  was  carried  out  successfully  by  means 
of  the  state-space  model  for  the  inverse  system.  For  the 
construction  of  the  LSV  model,  the  impulse  responses  of 
the  blurring  phenomena  were  approximated  in  degenerate  form 
via  series  expansions  using  orthogonal  functions.  For  details 
regarding  this,  see  [5.21].' 


5.53 


The  possibility  of  using  the  state-space  model  of  the 
inverse  system  for  the  restoration  of  the  original  object 
from  a  LSV  blurred  image  in  the  presence  of  additive  noise  is 
currently  under  investigation.  A  causal,  discrete  counterpart 
of  the  integral  equation  in  Phillips  [5.22]  was  implemented, 
and  for  variou'  signal- to-noise-ratios  (SNR)  the  deconvolution 
was  performed.  It  is  interesting  to  note  that  in  this  case, 
there  were  no  large  oscillations  and  as  the  SNR  increased  the 
difference  between  the  original  sequence  and  the  deconvolved 
one  decreased.  Also  the  deconvolution  of  motion  blurred  and 
coma  blurred  images  corrupted  by  noise  was  performed  for 
different  SNR.  The  restored  objects,  for  a  SNR=100  or  larger, 
were  easily  recognized  in  the  case  of  objects  with  well 
defined  edges,  such  as  a  white  letter  on  a  black  background. 
The  reconstruction  of  the  original  object  was  very  poor 
in  the  case  when  the  object  did  not  have  sharply  defined  edges 
Since  the  state- space  model  developed  works  very  effi¬ 
ciently  to  deblur  images  affected  by  2-D  linear  shift- 
varying  blurs,  its  use,  in  presence  of  noise  needs  to  be 
examined.  An  obvious  approach  would  be  to  filter  out  the 
additive  noise,  and,  subsequently,  obtain,  recursively,  the 
restored  object  using  the  state-space  model  of  the  inverse 
system,  already  developed.  Specifically,  let  x(k^,k2)  be 
the  blurred  image  with  additive  noise, 

x(k^,k2)  =  s(kj^,k2)  +  n(k^,k2)  (5.95) 

where  s(kj^,k2)  is  the  spatially-variant  blurred  image  and 
n(k^,k2)  is  the  additive  noise.  We  want  to  filter  out 
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n(k^,k2)  or  in  other  words  we  want  a  satisfactory  estimate, 
s(k^,k2),  of  the  blurred  image.  From  this  obtained  estimate, 
s(k^,k2),  the  deconvolution  can  be  implemented,  recursively, 
via  the  developed  state-space  model.  Some  assumptions  are 
necessary  before  n(k^,k2)  may  be  filtered  out  to  satisfaction. 

A  stationary  noise,  {n(k^,k2)},  which  is  uncorrelated  with 
{s(k^,k2)}  and  having  known  mean  as  well  as  correlation  may 
be  assumed.  The  sequence,  {s(k2^,k2)}  originates  from  a  spati¬ 
ally  blurred  object  and,  therefore,  it  is  inherently  non¬ 
stationary.  We  will  assume  that  {s(k^,k2)}  is  nonstationary 
in  the  mean  and  also  in  the  autocorrelation.  Under  these 
assumptions,  we  feel  that  a  solution  to  the  filtering 
problem  could  be  obtained  and  a  possible  approach  is  outlined. 

It  is  possible  to  transform  {x(k^,k2)}  into  an 
approximately  stationary  process,  {x(k^,k2)},  given  by  [5.23], 

x(k^,k2)  =  s(k^,k2)  +  n(k^,k2)  (5.96) 

where  {¥(kj^,k2)}  and  {n(kj^,k2)}  are  stationary  and  uncorrelated 

and  the  mean  and  correlation  of  {n(k^,k2)}  are  computable 

from  those  of  {n(kj^,k2)}.  From  the  process,  {x(k^,k2)}, 

an  estimate,  s(kj^,k2)  ,  of  s(k^,k2)  is  obtained.  This  is 
possible  to  do  via  use  of  Wiener  filtering  theory.  Next, 

an  inverse  transformation  to  that  employed  in  order  to  arrive 

at  (5.96)  is  applied  to  {J(k^,k2)}  and  an  estimate  s(k^,k2) 

of  the  original  image  s(kj^,k2)  is  obtained.  It  is  pointed 

out  that  if  the  noise  is  nonstationary,  further  assumptions 

are  necessary  (like  local  stationarity)  before  a  satisfactory 

solution  to  the  problem  is  expected. 
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Linear  Restoration  of  Bilinearly  Distorted  Image:  Several 
applications  require  restoration  of  bilinearly  (a  special 
type  of  nonlinear  map)  distorted  images.  Some  of  these 
applications  and  a  procediare  to  restore  bilinearly  distorted 
image  in  the  presence  of  additive  noise,  by  linear  filtration 
is  considered  in  [5.24],  It  is  also  known  how  1-D  bilinear 
transformation  (shift- invariant  or  special  shift-variant)  can 
be  computed  by  use  of  2-D  linear  optical  processors  [5,25]. 

It  has  also  been  pointed  out  in  [2,  p.  219]  that  properties 
of  n-D  bilinear  systems  can  be  inferred  from  the  investigation 
of  similar  properties  in  a  2n-D  linear  system.  These  inter¬ 
relations  between  a  bilinear  and  a  higher  dimensional  linear 
system  suggest  the  necessity  of  investigating  into  the 
possibility  of  restoring  bilinearly  distorted  images  using 
the  developed  state-space  model  (of  appropriate  spatial 
dimension)  for  linear  shift  variant  systems. 


5.56 
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g  4  [g(0),g(i).*  •  vg{N-i)]''‘  (3.i0) 

hen,  the  bilinear  system  description  in  (1.1),  or  (3.2),  can  be  written  as 

Bx  *  g  (3.11) 

herefore,  the  image  restoration  problem  of  interest  here  requires  the  finding  of  the 
onnegative  vector  f  which  satisfies  (3.7),  (3.8)  and  (3.11). 

Fact  3.1  :  Let  h^^^,  n=0,1,*  •  '.N-l,  B,  f,  and  g  be,  respectively,  as  in  (3.4),  (3.5),  (3.6), 
J.3),  and  (3.10)  and  h^^^  n=0,1,*  •  vN-l  and  f,  g  be  nonnegative.  Define  an  Nxl  vector  w, 
rhose  elements  are  given  by. 


(w)i  - 

nd  also  define  an  NxN  matrix  U 
■h{0)T  ' 

u  A  • 


(3.12) 


(3.13) 


h(N)T 


f  is  a  solution  of 


Uf  -  w 


(3.14) 


hen  X  obtained  from  f  via  (3.7)  and  (3.8)  will  be  a  solution  of 


Bx  -  g. 


(3.15) 


et  the  linear  subspace  R(B  )  be  defined  as 


R(B^)  4  {z  :  pB^y,  yeR*'^} 


(3.16) 


(/here  R^  denotes  the  N-dimensional  Euclidean  space.  Then,  the  projection,  x,  of  x  on 
l(B^)  will  satisfy 


<  (B)i^,x  >  -  (g)j,  i«1,2,'‘*,N 


(3.17) 


Proof  :  (3.14)  implies  that 
<h('~"'\f>  -  (w)j. 
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3.  Image  Restoration  From  Noncausal  Bilinear  Blurring 

In  this  section,  the  problem  of  restoring  an  image  which  is  blurred  by  the  bilinear 
system  described  by  (1.1),  will  be  discussed.  The  (NxN)  nonnegative  matrix,  is,  first 
defined. 


q(n;0,0) 

q(n;0,1) 

q(n;0,N-1) 

q(n;1.0) 

• 

q(n;1,1)  *** 

•  •  «  • 

q(n;1,N-1) 

• 

n-0,1,*  *  *,N-1 

(3.1) 

« 

q(n;N-1,0) 

• 

•  •  •  • 

q(n;N-1,1)  *  *  * 

• 

q(n;N-1,N-1) 

Then,  n-0,1, •  •  •,N-1,  are  symmetric  matrices  if  Y(mj,m2)  in  (1.2)  is  symmetric. 

Assume  that  Y(mj,mj)  is  symmetric  (4,6,7,  13l  throughout  this  section.  With  in 


(3  1),  aquation  (1.1)  can  be  rewritten  by 

g(n)  -  f^Q<"V  n-0,1, ••  vN-l  (3.2) 

where 

i  -  (f(U),f(1),***,f(N-1)f  (3.3) 

If  the  system  is  completely  coherent,  then 

Q{n)  ,  n-0,1,*  *  VN-I  (3.4) 

where 

h('’)  -  [h(n;0),h(n;1),*  •  •,h(n;N-1 )]''’,  n-0,1,*  •  •,N-1  (3.5) 

Let  (A)j  denote  the  i^^-row  of  A.  Let, 

(B),  A  [(Q('"1))i,(Q('"‘'))j,***,(Q(*‘‘'))pj],  i-1,2,***,N  (3.6) 

X  -  ff^  (3.7) 


and 


X  -  [(X)j,(X),.***.(X),^f  (3.8) 

Then,  using  the  usual  Euclidean  innerproduct  notation,  <*,*>,  (3.2)  is  expresseible  as, 

g(n)  -  <  (B)p,,^^^,  X  >,  n-0,1,*  * ‘.N-l  (3.9) 


Let 
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has  been  done  in  [17],  For  a  (31x31)  image,  Rg  <  0.0065  guarantees  the  occurrence  of 
the  desired  effect  and  for  each  segmented  region,  the  support  of  hj.(n^,n  j;mj,mj)  is 
obtained  from: 


hQ(n  ^  #n  2f  ni  ^  ,m  2} 


nonzero  ,  m j  <  n ^  and  m2  ^  n 2 
0  ,  otherwise 


(2.21) 


Then,  each  segmented  image  is  blurred  bilinearly  by  using  (2.13).  Figure-5(a)  shows  the 
original  31x31  image  and  Figure-5(b)  shows  the  blurred  image.  To  simulate  the  blurred 
image,  A  in  (2.19)  has  been  chosen  as  0.254.  Also  y  2(^0  z-^  1-^2)  9'ven  by 

Y2(m^.m2.4i,2.2)  *  sinc((mj-il J/Nf,]sinc[(m2-il2)/Nf,]  (2.22) 

where  Nf,  is  the  maximum  possible  support  width  of  hj.(nj,n2;mj,m2)  for  any  {nj,n2). 
The  system's  coherent  impulse  response  is  given  by 


h2(n  j.n2;m2.m2)  ■  /H^Irrpr^lrnTrnTT  (2.23) 

Each  function  has  been  chosen  to  ensure  the  nonnegativity  of  the  resulting  OIR.  By 
applying  the  algorithm  in  Figure- 1  to  each  segmented  output  quadrant  image  and  by 
combining  those  results,  the  original  image  is  recovered  as  in  Figure-5(c). 
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The  recursive  implementation  analogous  to  (2.8),  (2.9)  and  (2.W)  also  can  be  obtained. 
Since  (2.13)  characterizes  a  quarterplane  2-D  discrete  system,  row-by-row  recursion,  i.e., 
(Oj.nj)  (nj  +  l.nj)  •••  -*•  (N~l.nj)  -►  (O.n^  +  I)  (Tn^  +  I)  -►  •  •  •,  and  the 

column-by-column  recursion  are  two  of  several  possibilities  for  implementing  the 
recursion.  The  next  equation  provides  the  input/output  description  of  2-0  discrete 
bilinear  system  whose  DIR  has  support  in  a  nonsymmetric  half-plane  on  (m^,m2)-piane 
or  (a^.J-jj-plane. 


V(nj,n2) 


"i  N-1  "i  N-1 

^  ^q2(nj.n2;m^,m2,JLi,i2)u(mj,m2)u()lj,)l2) 

m ,  *0m  ,=02, ,  ■0£,»0  *  *  2  i  z  i  z  1  z  i  z 


(2.17) 


or. 


J-1  "2 


S  2 

mj*0m2*02,^*02,2*0 


N-1  "2 

£  Z  .q2(nj,n2;mi.m2,2,i,2,2)u(mj,m2)u(2,j,2,2)  (2.18) 


Yini.n^) 

(2.17)  can  be  solved  recursively  column-by-column,  and  (2.18)  row-by-row 


To  simulate  the  proposed  algorithm,  we  synthesize  a  special  type  of  noncausal  bilinear 
system  by  employing  a  coma  type  lens  aberration  as  the  system’s  coherent  impulse 
response.  The  coma  aberration  is  described  in  [17,  23].  It  will  be  briefly  described  here. 

In  the  rectangular  coordinate  system,  it  is  easy  to  explain  this  coma  aberration  by 
segmenting  the  input  support  region  into  4  quadrants  as  in  Figure-3.  And  for  each 
segmented  quadrant  the  index  will  be  reordered  the  center  point  of  the  original  input  as 
the  origin,  i.e.,  (0,0),  of  the  each  segmented  image.  In  Figure-3,  <m^,m2>  denotes  the 
original  index  of  (31x31)  original  image  and  (m^,m2)  the  new  index  for  each  segmented 
part.  The  point  spread  function,  hQ(n  ^,n2;m  ^,m2),  of  the  lens  coma  aberration  is  given 
by  [17] 

h^.(n^,n2,■m^,m2)  -  [1/{(m ^ A)*+(m2A)^}]h,(x^,X2)  (2.19) 

where  A  denotes  the  sampling  distance  and 

Xj  4  (m^n^+m2n2)/(m^  *+m2*) 

X2  A  (mjn2-m2nj)/(mj  *+m2*) 

ho(Xi,X2)  -  |(2C)/y^7^^3v*’  .  (x^,X2)cl  (2.20) 

(  C//Xj*-3x2*  ,  (x^,X2)ell 

and  regions  I  and  II  are  defined  in  Figure-4.  From  Figure-4,  it  is  clear  that  by 
constraining  to  be  less  than  a  certain  value,  it  is  possible  to  ensure  that  the 
i^^-quadrant  input,  i*1,2,3,4,  does  not  affect  the  j**^-quadrant  (i¥j)  output.  This  analysis 


Solve 


AnCf(n)l"  ♦  Bn[f(n)]  - 


Figur8-2. 


Added  Algorittim  for  Real  DIR 


•viV 
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Figure-1.  Recursive  Restoration  Algorithm 
for  Nonnegative  DIR  Case 
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Suppose  the  values  for  f(m),  m=0,1,*  • ‘.n-l,  are  known;  then,  Ap,  Bp  and  Cp  can  be 
computed  by  using  (2.5),  (2.6)  and  (2.7).  Hence,  f(n)  may  be  obtained  by  solving  (2.12). 
Since  all  the  coefficients,  Ap,  Bp  and  Cp,  in  (2.12)  are  nonnegative,  the  equation  (2.12)  has 
two  real  solutions,  one  of  which  is  nonnegative  and  the  other  is  negative.  The 
nonnegative  solution  is  assigned  to  f(n)  and  the  recursion  is  continued.  The  flowchart  for 
this  is  given  in  the  Figure-1.  In  the  Figure-1,  the  initialization  S,=-1  will  not  be  necessary 
in  this  case.  After  recursion  is  completed,  the  desired  nonnegative  output,  f(m), 
m=»0,1,»  •  •,N-1.  is  uniquely  obtained. 


Next,  suppose  that  the  elements  of  the  DIR,  are  not  all  nonnegative.  Then,  the 
inequalities  (2.3)  and  (2.4)  may  not  hold.  To  proceed,  it  is  assumed  that  q(0;0,0)>0.  The 
flowchart  in  Figure-2,  when  embedded  appropriately  in  the  hatched  box  of  Figure-1 
provides  a  brute-force  implementation  of  the  solution  f(m),  m=0,1,*  •  •,N-1,  when  this 
soiution  exists. 

For  a  2-0  discrete  bilinear  system,  the  counterpart  of  the  expression  in  (2.1)  is 

n^njn^nj 

Y(n,,n2)  “  2  S  2  2  q2(hj<n2;m (2.13) 

‘  *  mi=0m2»01lj=0fi,2»0  z  ‘  i  i  i  ‘  i  i 

where  all  the  notations  are  self-evident.  The  counterparts  of  (2.5),  (2.6),  (2.7)  in  this  case 

are  given  below. 


^  q 2(h i/n2,n j,n j.n j,n j) 

n  j-1  Oj-l 

Bn  n  4  2  „{q2(^.,n,;m,,m2,n.,^,)+q,(n^,n2;n^,n2,m^,m2)}u(m^,m^) 

*  mj-0m2»0 


(2.14) 


•^2 

m- 

i'" 

2;m,n 

2,n  j,n 

"2 

-1 

+  2 

m- 

;,n 

j/H  j,m,n 

"i- 

-1  02- 

in,- 

,0  2 

4  v(n 

i,n 

2-2 

m^ 

2 

-0m2  = 

■0^1 

"1 

-1  02- 

1  n 

1-1 

-2 

m, 

.-Om.- 

'o£ 

(Hi.Oj 

;;m  j,i 

(2.15) 

a,) 


n^-l  nj-l  nj-1 


^2  ^{q,(n.,n,;m.,m..n.,Jl)+q2(n.,n2;n.,Jl,m.,m2)}u(m.,m,)u(nj,a) 
mj«0m2“04*0 

n^-ln^-l  nj-l 

^  q,(n  ,.n,;m,n,,l,n2)u(l.n,)+2  q 2(0  , ,n  ,;m,n ,,n  .  ,l)u(n  .  ,i)]u(m,n  j) 

m-0  4»0  z  ‘  z  z  i  i  t  I  *  t  1 


m 

n  2“1n  j -1 


02-1 


-2  [2  q,(n  ,,n,;n  ,,m,i,n,)u(a,n,)+2  ^q  ,(n  .  ,n  2;n  .  .m.n  .  ,2,)u(n  ,2,)]u(n  .  ,m)(2.16) 

m»0  1-0  ^  121  z  z'  z  1  z  1  1  I  I 
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2.  Image  Restoration  From  Causal  Bilinear  Blurring  (Noise-Free  Case) 
An  1-0  discrete  causal  bilinear  system  will  be  represented  by. 


z  2 


q{n;m  ^,m2)f(m  Jf(m^). 


The  case  of  nonnegative  OIR  will  be  considered  first.  Suppose  that  we  are  given  g(n), 
n=0,1,*  • ‘.N-l  and  q{n;mj,m2),  0  <  m^^.m^  <  n,  n*0,1,*  •  vN-l,  and  these  are 
nonnegative.  The  nonnegative  input  f(m),  m»0,1,«  • '.N-l,  is  to  be  recovered.  This  can  be 
done  by  rewriting  (2.1)  as 

g(n)  =  q(n;n,n)f^(n)  +  (2  ^^{q(n;m,n)+q(n;n,m)}f(m)]f(n) 

+  22  q(n;mj,m2)f(mj)f(m2).  (2.2) 

m^=0m2=0  i  i  I  i 

Since,  for  the  noise-free  case,  q(n;m^,m2)  >  0  and  f(m)  >  0,  therefore 

q(n;n,n)f ^(n)  +  [2  ^  {q(n;m,n)+q(n;n,m)}f(m)]f(n)  >  0.  (2.3) 

m»0 

From  (2.2)  and  (2.3),  the  inequality  in  (2.4)  follows. 

g(n)  -  2  \2  ^  q(n;m.,m2)f{mi)f(m2)  >  0.  (2.4) 

mj=0m2=0  ‘  i  — 


Af,  4  q(n;n,n) 

Bp  4  2  ^  {q(n;m,n)+q(n;n,m)}f(m) 


Cn  A  =  g(n)  ■  2  ^  q(n;mj,m2)f(mj)f(m  ) 

"  "  mj-0m2=0  I  4  1  i 

Cjn)  can  be  obtained  by  implementing  the  following  recursion. 


cp  -  0.  k-0,1,*«vn 

c|^^  *  +  (^  ^  {q(n;m,k-1)+q(n;k-1,m)}f(m)]f(k-1) 

+  q(n;k-1,k-1)f^(k-1),  k*1,2.***,n 
Cj">  -  g(n)  -  Ci,"> 

From  (2.4),  (2.5),  (2.6)  and  (2.7) 


(2.10) 


An  ^0-  Bn  ^  0'  Cp  >  0. 

The  expression  (2.2)  can  be  rewritten  in  terms  of  B^,,  and  Cp,  by 
A^f^n)  +  Bnf(n)  -  C^  -  0. 


(2.11) 


(2.12) 
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h(n;m)  «  nonzero  ,  m  <  n 

.  0  ,  otherwise, 

then  it  is  easy  to  see  that  the  system  OIR  will  be  of  the  form 


(1.12) 


q{n;m^.m2)  »  j  nonzero  ,  m^  <  n  and  m^  <  n  (113) 

I  0  ,  otherwise. 

For  the  incoherent  imaging  system  of  patterns  on  the  translucent  scattering  substrates, 
the  causality  of  the  lens  irradiance  spread  function  and  the  substrate  scattering  function 
will  result  in  the  system  DIR  of  the  form  in  (1.13).  A  physical  realization  of  the  causal 
scattering  substrate  has  not  yet  been  reported.  However,  the  example  of  the  lens 
irradiance  spread  function  of  the  type  in  (1.12)  can  be  formed  as  we  shall  see  in  section  2 
by  considering  the  lens  coma  aberration  [17],  Since  the  irradiance  spread  function, 
hj.(n;m),  of  the  lens  is  related  to  the  system's  coherent  impulse  response,  from(1.4),  by 

hjn:m)  =  |h(n;m)|^,  (1.14) 

one  can  say  that  the  support  of  h(n;m)  is  the  same  as  that  of  hj,(n;m). 

In  this  paper,  we  will  mainly  study  the  original  image  intensity  restoration  from  the 
image  intensity  blurred  by  the  system  described  by  (1.1).  For  the  restoration  of  the 
bilinearly  blurred  images,  nonlinear  [18,  19]  as  well  as  linear  methods  [12,  20,  21,  22]  have 
been  considered.  The  linear  methods  used  for  the  restoration  of  a  bilinearly  blurred 
image  are  successful  in  the  incoherent  case,  but  the  restoration  becomes 
poor  [18,  20,  21,  22]  as  the  blurring  system  approaches  the  completely  coherent  case. 
The  nonlinear  method  [19]  has  been  successfully  applied  to  restore  only  those  images 
which  are  of  low  contrast  following  blurring  by  a  coherent  system. 


A  procedure  to  restore  images  blurred  by  completely  coherent  systems  is  described  in 
section  3.  In  several  applications  including  the  imaging  system  of  patterns  on  the 
translucent  scattering  substrate  [13]  and  a  poaion  of  the  diffraction-limited  imaging 
system  [7L  the  characterizing  OIR  is  nonnegative.  Except  in  section  2  this  nonnegativity 
property  of  the  OIR  is  assumed  to  hold.  In  section  2,  a  recursive  scheme  for  restoring  a 
class  of  bilinearly  blurred  images  is  developed.  In  section  3,  a  technique  for  restoring  a 
wider  class  of  bilinearly  blurred  images  is  considered.  In  section  4,  the  effect  on  the 
quality  of  restoration  due  to  the  presence  of  additive  noise  is  analysed. 
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g(n)  “  Z  |h(n;m)|^lf(m)|^  (1.4) 

m*0 

where  g(n)  and  |f(m)|^  represent,  respectively,  the  image  intensity  at  n  and  the  object 
intensity  at  m.  For  the  completely  coherent  case,  in  which  the  field's  coherence  function 
is  given  by  [4,  6] 


^(mj.mj)  *  1  .  for  all  and  m2  (1.5) 

we  have 

g(n)  -  jS’\(n:m)f(m)|*.  (1.6) 

m-0 

The  region  falling  between  the  two  extremes  of  complete  incoherence  and  complete 
coherence  is  the  region  of  partial  coherence  and  the  partially  coherent  system  is 
described  by  (M)  and  (1.2).  When  the  following  condition  holds, 

yCm^.mj)  -  Y*(»n2.mJ  (1.7) 

the  field's  coherence  function  is  said  to  be  symmetric  and  in  that  case  the  OIR 
q(n;mj^,m2)  is  also  symmetric.  This  implies  that, 

q(n;mj.m2)  •  q*(n;m2.mj).  (1.8) 

When  (1.8)  holds  and  the  input  sequence  is  real,  (1.1)  is  a  hermitian  form  and,  therefore, 
can  be  written  as 

g(n)  -  S  ^  Relq(n;m.,m2)lf(m,)f(m2).  (1.9) 

m^*0m2*0  *  *  * 

Hence,  the  imaginary  part  of  the  system's  DIR  doesn't  contribute  to  the  image  formation. 

In  the  imaging  system  of  patterns  on  translucent  scattering  substrates  [13],  the  OIR  can 
be  given  by 

q(n;mi.m2)  ■  (1/2l[M(mj,m2)S(n;mj)+M(m2,mj)S(n;m2)]  (1.10) 

where  M(m^,m2}  represents  the  substrate  scattering  function  at  surface  point  m^ 
resulting  from  a  small  illumination  spot  of  unit  power  at  surface  point  mj  and  S(n;m)  is 
the  lens  irradiance  spread  function.  From  (1.10),  we  have 

q(n;mj,mj)  »  q(n;m2,mj)  (1.11) 

and  q(n;m^,m2)  always  takes  the  real,  positive  value  [13]. 


For  the  partially  coherent  imaging  system,  if  the  system's  coherent  impulse  response 
has  the  following  support 


1.  Introduction 


A  1-0  continuous  bilinear  transformation  arises  in  the  second  order  term  of  the  Voiterra 
series  representations  of  the  nonlinear  system  [1,  2,  3,  4].  In  the  discrete  domain,  the  1-0 
bilinear  transformation  with  the  finite  support  range  input  can  be  described  as 

g(n)  -  S  ^  E  ^  q(n;mi,m2)f(m^)f(m2)  (1.1) 

mj=0m2“0  i  ^  1  i 

where  {g(n)}  and  {f{m)}  are,  respectively,  the  output  and  input  sequences,  assumed  real, 
q(n;mj,m2)  represents  the  system  response  at  the  output  coordinate  n  due  to  two 
impulses  at  the  input  coordinates  m=mj  and  m=m2,  and  N  is  the  size  of  the  input 
support. 

Recently,  in  the  optical  image  processing  area,  this  bilinear  transformation  has  been 
studied  extensively  to  analyze  the  optical  imaging  system.  These  areas  include  the 
partially  coherent  imaging  system  [4,  5,  6,  7,  8,  91  the  magnification  type  X-ray  imaging 
system  [10,  11,  12],  and  the  incoherent  imaging  patterns  on  translucent  scattering 
substrates  [13].  Ref[4]  describes  the  general  properties  of  the  optical  bilinear 
transformation  in  detail.  The  relationship  between  the  1-D  bilinear  transformation  and  the 
2-D  linear  transformation  has  been  studied  in  [4,  14,  15,  16]. 

In  the  partially  coherent  imaging  system  and  the  projection  type  imaging 
system  [10,  11],  the  double  impulse  response  (DIR),  q(n;mj,m2),  of  equation  (1.1)  can  be 
represented  by 

q(n;mj,m2)  -  h*(n;mj)h(n;m2)Y(mj,m2)  (1.2) 

where  h  (n;m)  denotes  the  complex  conjugate  of  h(n;m),  h(n;m)  represents  the  system's 
coherent  impulse  response  at  n  due  to  the  impulse  at  m  and  Y('Tii,m2)  is  the  field's 
coherence  function.  The  term  Y(fT’i'nij)  represents  the  correlation-like  coefficient 
between  the  object  intensities  at  m^  and  m2.  The  derivation  of  the  expression  (1.2)  and 
the  further  detaiied  informations  about  Y(ttij^,m2)  can  be  obtained  in  [6,  7,  8]  and  various 
other  related  materials. 

For  the  completely  incoherent  case,  the  light  from  each  point  in  the  object  is  assumed 
to  be  statistically  independent  of  light  from  every  other  point  [6].  Therefore,  in  this 
case  [4,  6], 


Y(m^,m2) 


1  ,  m,*m, 

0  ,  otherwise 


(1.3) 


and,  from  (1.1), (1.2)  and  (1.3),  we  obtain 
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BILINEARLY  BLURRED  IMAGE  RESTORATION 


l« 


On  squaring  both  sides  of  the  preceding  equation,  and  then  using  successively  (3.12),  (3.4) 
and  (3.11),  one  obtains 

(g)i  -  f^h('"1)[h('“‘')rf 

-  fTQ(i-1)f 

-  <(B)j''',x>  (3.18) 

Therefore,  if  f  is  a  solution  of  (3.14),  then  g  is  a  solution  of  (3.15). 

To  prove  (3.17).  x  is  decomposed  as, 

X  =  X  +  X  (3.19) 

where  xeR(B^)  and  xeN(B),  where  N(B)  is  defined  as, 

N(B)  4  {y  :  By»a  (3.20) 

since  R(B^)  in  (3.16)  and  N(B)  ar«  orthogonal  complementary  spaces  with  respect  to  R*'^*. 
Premultiplying  B  on  both  sides  of  (3.19), 

Bx  =  Bx  +  Bx  •  Bx  (3.21) 

But,  by  (3.15)  the  left-hand  side  isg.  Hence, 

Bx  -  g  (3.22) 

and  X  satisfies  (3.17). 

The  above  Fact  says  that  if  the  each  row  of  the  matrix  B  satisfies  (3.4),  (3.5)  and  (3.6), 

then  by  forming  a  nonnegative  vector  w  from  g  as  in  (3.12)  and  solving  (3.14)  for  f,  one 
can  restore  the  original  image.  When  the  nonnegative  vector  w  is  formed  by  (3.12),  there 
are  2*^^  possible  ways  to  sat  the  value  (w)j  from  (g),-,  i=1,2,  •  •  vN.  But  the  negative  sign  of 
any  elements  in  w  will  violate  the  necessary  condition  for  f  to  be  nonnegative  in  (3.14). 
Hence,  w  may  be  formed  uniquely  from  g  without  violating  the  necessary  condition  that  f 
be  nonnegative,  by  using  (3.12).  If  h^'\  i“0,1,*  • ‘.N-l,  are  linearly  independent,  the 
nonnegative  solution  f  will  be  unique.  For  the  completely  coherent  case,  the  matrix  B  in 
(3.11)  satisfies  the  condition  in  the  Fact  3.1  by  (1.2),  (1.5)  and  (3.1).  Therefore,  one  can 
exactly  recover  the  original  image  for  this  case  provided  that  h^’\  i=0,1,*  •  •,N-1,  are 
linearly  independent. 

^  T 

Due  to  the  Fact  3.1,  if  the  projection,  x,  on  R(B')  of  x  is  known,  where  B  can  be 
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represented  by  (3.4),  (3.5)  and  (3.6),  then  the  null  space  component,  x  of  x  can  be  obtained 
from  X  so  that  x  can  satisfy  (3.7)  and  (3.8).  However,  for  the  partially  coherent  system, 
the  matrix  B  generally  does  not  satisfy  the  conditions  in  (3.4),  (3.5)  and  (3.6).  The  Fact  3.1 
will  be  adapted  for  use  in  the  partially  coherent  case,  by  approximating  the  basis  of  R(B^) 
to  the  form  in  (3.4),  (3.5)  and  (3.6).  To  do  this  a  few  terminologies  and  well-known 
Lemmas  will  be  introduced. 


Definition  3.1  :  The  Euclidean  matrix  norm  of  the  matrix  D  will  be  denoted  by  |ID||  and 
defined  as 

IlDll  4  /TrlD^O]  (3.23) 

where  TrA  denotes  the  trace  of  the  square  matrix  A. 


Definition  3.2  ;  Let  S  be  the  set  of  the  matrices  satisfying 


S  »  {D;  D-dd'*’,  deR'^} 


(3.24) 


!• 


Then,  we  mean  by  'the  best  approximation  of  a  matrix  Q  on  S'  the  matrix  D*  which 
satisfies 


llQ-ol 


min  ||Q-D|I 
DeS" 


(3.25) 


The  above  two  definitions  can  be  explained  in  other  way.  Let 

a  >  I(Q)j,(Q)2,««-,(Q)^,f  (3.26) 

e  -  [(D)j,(D)2,«--,(D),^jf  (3.27) 

where  D*dd^,  deR*'^.  Then,  the  best  approximation  of  a  matrix  Q  on  S  is  seeking  the 

it 

vector  e  which  has  the  form  in  (3.27)  and  minimizing  the  usual  Euclidean  vector  norm  of 
[a-e]. 


Definition  3.3  :  A  matrix  E  is  said  to  be  idempotent  if  E^“E. 

The  following  Lemma  3.1  is  well-known  [26]  and  its  proof  will  be  omitted. 

Lemma  3.1  :  Let  Q  be  a  NxN  symmetric  matrix  with  distinct  real  eigenvalues, 
Then,  the  idempotents  Ej,  i*1,2,*  •  •  ,t  having  the  properties 

1)  EjEj  - 


0  if  i¥j 
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I 

’  'n 

3)  Q  = 

exist,  where  Ijj  denotes  the  NxN  identity  matrix.  Moreover,  Ej,  i=1,2, •  •  •,r  are  unique  and 
Ei=yi\/i^,  where  Vj  is  the  eigenvector  of  Q  corresponding  to  Xj,  i=1,2,»  •  vr. 

Lemma  3.2  :  If  E  is  idempotent,  rank  E  =•  Tr  E. 

Proof  :  Since  E  is  invariant  under  the  square  operation,  the  eigenvalues  of  E  are  1  and  0. 

Hence,  the  multiplicity  of  the  1  is  the  rank  of  E.  But  the  sum  of  the  eigenvalues  of  any 
matrix  is  the  trace  of  the  matrix.  Therefore,  Tr  E  =  rank  E. 


I* 


Then,  we  state  the  Theorem  3.1  and  present  its  proof. 

Theorem  3.1  :  Let  S  be  the  set  as  in  (3.24).  Then,  the  best  approximation,  D*  of  a 
symmetric  nonnegative  NxN  matrix  Q  on  S  is  given  by 

D*  =■  X^Ej  =*  (3.28) 

where  X^  and  y^  are,  respectively,  the  dominant  eigenvalue  and  the  dominant 
eigenvector  of  the  matrix  Q  and  E^  denotes  the  dominant  eigenspace  of  Q. 

Proof  :  Suppose  the  rank  of  Q  be  r.  Let  the  positive  eigenvalues  of  Q  be  ordered  in 
such  a  way  that 

Xi>  X2>  •••>  Xp  (3.29) 

where  p  is  the  number  of  the  positive  eigenvalue,  and  the  negative  eigenvalues  in  such  a 
way  that 

Vl,^  ^p+2^  ^r  (3-30) 

Then, 


. 


Xr+i»*  • '-Xivj-O  (3.31) 

Now,  suppose  that  yj  denote  the  orthonormalized  eigenvectors  of  Q  corresponding  to  X,- 
i=1,2, •••,N.  Then.  yi.***.y|\|  form  an  orthonormal  basis  of  Let  d  be  the  vector 

which  minimizes 

llQ-dd'’’!!  (3.32) 

Then,  d  can  be  expressed  as  a  linear  combination  of  vj,  i  =  l,2, •  •  vN. 


Since 


-  TrQ^  -  Tr[Q(dd''‘)l  -  Tr((dd‘'')Q]  +  Tr( 


Tr[Q(dd'’’)]  =  Tr[(dd''')QL 

|lQ-dd''’||*  =■  TrQ^  -  2Tr[Q(dd''')]  +  Tr(dd''‘)^ 
By  the  Lemma  3.1, 

TrQ^  *  Tr(Z  =  Trfz 

i=1  ‘  i=1  '  ' 

Hence,  by  the  Lemma  3.2, 

TrQ^  =  Z  Xi^TrE:  =  Z  Xj^ 
i-1  '  '  i=1  ' 

For  the  second  term  in  (3.34), 

TrCOyd'^)]  . 

-  Tr(^_^Xia|yi)(|_^oi|,vJ) 

Since 


TrlViv/l 


bi-Ykl 


1,  i=k, 

0,  otherwise. 


Tr[Q(dd‘’’)]  =■  ^  Xidj^ 
i=1  '  ' 

For  the  third  term  in  (3.34), 

Tr(ddT)k  .  Trt(?,0|yi)(?^aiy|T)<2_^Okyk)(2_^a|iyiT„ 


^  ^  £  ajaj^aaTrlvjVo^] 

i=1j=U=>1  '  J  ^ 
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By  substituting  (3.35),  (3.36)  and  (3.37)  into  (3.34),  one  obtains 

|lQ-ddT|l!  .  V(Xi-a|")*  «  [!  Xi'  .  ?  J  ,ei|‘dk‘(i-«|K)I 

1*1  '  '  i*r+l  '  i»1k*1  J 


(3.37) 


A  £,(V«i*)*  ^  0 


N  N  N 

t,  A  Z  a/*  +  E  t  >  0 

^  *■  i=rf1  '  »1k»1  i  ^  - 


Then, 


llQ-dd^ll" 


e  ^  can  be  rewritten  by 


=  S  (Xj-aj*)*  +  Z 


(3.38) 


(3.39) 


(3.40) 


(3.41) 


a  *  [tti  >  0 


(3.42) 


By  varying  a  from  [0,0,  ••  *,0]^,  one  can  sot  the  a*  which  minimizes  Ej+e^.  Starting 
from  [0,0,  •  •  •  ,0]^,  E^  is  nondecreasing  function  of  g  as  ll^i  increases.  However,  Ej  is 
decreasing,  monotonically  as  ||a||  increases  from  [0,0,  •  •  vOl^  until  a  reaches 

a  »  [X^,X2,*  •  •,Xp,0,*  ••,0]^ 

Then,  e  ^  is  increasing  monotonically  as  ||gj|  increases.  Based  on  these  facts,  one  can 
conclude  that  the  optimal  g  minimizing  ||Q-dd  |)  exists  within  the  region  formed  by 
[0,0,*  •  *,0]^  and  [X  i,X2.*  •  •,Xp,0,0,*  •  vO)^. 

Suppose  g<^^X^,0,*  01^.  Then, 

*  5:  Xj2-x  ^ 

It  i.l  •  i 


-  [o,x,,o,* •  *,0]^. 


Then, 


Cd,‘s:,lW  .  J,Xi‘-X/. 

By  checking  in  this  way,  one  obtains 
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Hence,  so  far.  is  attained  as  the  minimum.  Now,  let 

Since 


.  Then, 


Xj*  <  XjXj. 


^(k,i)  «  [0.0,*  • ‘.OA  1^.0, •• ‘.O.X 5^,0,*  •  •.oF,  1<  lc,£<  p 

For  the  similar  argument  can  be  applied.  By  continuing  these  procedures  until 

o  •  IXj.Xj.-'vXp.O.-'vO]'*’ 

IS  reached,  one  can  show  that 


■  min({Q-dd^(|*. 

Therefore, 


Oj  *  /X7,  Oj-a,**  •  •»apj*0  (3.43) 

will  give  the  optimal  vector  d  minimizing  l|Q-dd^l|.  But, 

dd^  -  (/rjyj)(/X'iyj^)  *  ^lEj. 

Hence,  the  Theorem  3.1  has  been  proved. 


The  error  resulted  from  the  above  approximation  is  given  by 

|lQ-dd‘'‘lI  •  jZ^Xj^  (3.44) 

By  applying  the  Theorem  3.1  to  i='0,1,»»«  ,N-1  in  (3.1),  one  will  get  i*0,1,*** 
,N-1  which  is  nonnegative  due  to  the  Perron-Frobenius  Theorem  [24].  That  is,  the 

A 

nonnegativity  of  the  system  matrix  B  which  is  the  approximate  version  of  B  in  (3.11)  is 

A 

maintained,  while  the  each  row  of  the  matrix  B  satisfies  the  conditions  (3.4),  (3.5)  and 
(3.6)  in  the  Fact  3.1. 


By  combining  the  Fact  3.1  and  the  Theorem  3.1,  the  image  restoration  algorithm  from 
the  partially  coherent  system  can  be  obtained  as  follows: 


Step  1)  Obtain  Nx1  vectors  d:,  i»0,1,*  •  ‘.N-l,  by 


^i  *  (3-45) 

where  Xjj  and  Vj^  denote,  respectively,  the  dominant  eigenvalue  and  eigenvector  of 
Form  N^xl  vectors  e-,,  i“0,1,*  •  vN-l,  by 


e,  -  [(□(')) 


0^'^  = 


(3.46) 

(3.47) 


Step  2)  Obtain  the  minimum  norm  solution,  x,  of  (3.11). 


Step  3)  Obtain  Nx1  vector  w  by 


1/<ei_^,x>  ,  >  0-  i*1,2,»**,N. 

0  ,  otherwise 


(3.48) 


Step  4)  Solve  NxN  linear  system  Uf**w  for  f,  where 


(3.49) 


I  J  ■ 

* 

Step  5)  If  the  solution  vector  f  has  negative  element,  make  it  zero.  Denote  it  f  .  Then, 

Hr 

f  is  the  approximate  solution. 

By  applying  the  above  algorithm,  one  wilt  have  the  squarred  error  of 


(■O  1*0  '  ' 

where  Xjj  is  the  eigenvalue  of  and  Vjj  is  the  eigenvector  of 


(3.50) 


This  error  results  from  approximation  of  the  range  space  of  B  by  a  space  satisfying 
(3.4),  (3.5)  and  (3.6).  Figure-6  will  help  understanding  their  relationships.  In  Figure-6,  the 
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matrix  E  denotes  the  matrix  resulting  from  the  approximation  of  the  matrix  B  in  (3.11)  by 
applying  the  Theorem  3.1.  The  dotted  curve  represents  the  nonlinear  space  which 
consists  of  the  vector  having  the  form  defined  in  (3.7)  and  (3.8).  Then,  solving  Uf=w  in 
the  step  4)  will  correspond  to  seek  a  vector  x  on  the  dotted  curve  so  that  the  projection 
of  X  on  R(E^)  be  OB.  But,  the  true  solution  00  is  on  the  intersection  of  AD,  which 
represents  the  set  of  all  (least  squares)  solution  vectors  of  (3.11),  and  the  dotted  curve.  If 
the  true  projection  of  OD  on  R(E^),  that  is,  OC,  is  given,  the  true  solution  OD  can  easily 
be  obtained.  For  the  completely  coherent  case  or  the  completely  incoherent  case,  R(B^) 
coincides  with  R(E^).  Hence,  the  true  projection  OC  can  be  obtained  only  by  computing 
the  minimum  norm  solution. 

The  innerproduct  operation  in  the  step  3)  includes  the  projection  operation  of  OA  on 
R(E^),  that  is,  obtaining  OB.  This  step  results  in  the  error  because  the  true  projection  of 
OD  on  R{e'*')  is  OC.  The  reason  why  the  innerproduct  operation  in  step  3)  includes  the 
projection  operation  is  following. 

The  projection  of  OA  on  R(E^)  is  represented  by 

OB  =  (E'*‘E)(0A) 

since  E'^E  is  idempotent  and  the  projection  operator  on  R(e'’')  [25].  Suppose,  now  we 
want  to  compute  the  value  (E)(OB)  to  establish  the  system  equation  in  the  form  of  Ex»g 
which  is  the  approximated  version  of  Bx=g.  Then, 

(E)(OB)  «  (E)(E'^E)(OA)  -  (E)(OA) 

since  EE'^E=E  [25].  Since  (E)(OA)  represents  the  innerproduct  operation  in  step  3),  this 
step  includes  the  orthogonal  projection  operation  on  R(E^). 

One  more  thing  to  be  noted  is  the  possible  negative  values  resulting  from  the 
innerproduct  operations  in  (3.48).  In  the  proof  of  the  Fact  3.1,  the  innerproduct  result  of 
(3.18)  never  becomes  negative  because 

<  (B)j,  X  >  -  <  (B)j,  X  >  -  [<h^'\d>]^  >  0 

But,  for  the  partially  coherent  case,  since  R{B^)  does  not  exactly  coincide  with  R(E^),  the 
first  equality  in  the  above  expression  doesn't  hold.  That  is, 

<  (E)j,  X  >  V  <  (E)},  X  >  +  <  (E)j,  X  > 
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since  x  may  contain  the  connponent  belonging  to  R(E^).  Hence,  if  the  error  in  (3.43)  is 
sufficiently  small,  then  x  can  be  said  to  belong  to  N(E).  Hence,  the  innerproduct  values  in 
step  3)  have  to  be  nonnegative. 


The  2-D  bilinear  system  representation  corresponding  to  (1.1)  is  given  by 


N-1  N-1  N-1  N-1 

Y(n  i-n  2)  -  2  ^  ,q2  i'"  2'*"  2'^  2)u(m  i<m  2)0(2,  ^ ,Z^). 

Initially  the  following  form  of  the  2-0  OIR  is  considered. 


(3.51) 


q2(nj,n2;mj,m2,2,^,2,2)  =  q(n  j;m  ^,2, 2)q(n2;m  2.2,2).  (3.52) 

In  the  Kohler  illumination  system  [7]  with  an  incoherent  square  source,  the  2-D  field's 
coherence  function  is  of  the  form. 


Y2{mi,m2,2,i,2,2)  ’  Y(fni,Jli)Y(m2,il2)  (3.53) 

where  Y(fT’,^)  “  Cj  •sinclC2  *(m-2,)]  and  and  C2  are  constants.  It  is  well-known  [7] 
that  the  Fraunhofer  diffraction  pattern  for  a  square  aperture  results  in  a  product  separable 
form  of  the  coherent  impulse  response,  i.e.,  the  2-D  counterpart  of  h(n;m)  in  (1.2)  is  of 
the  form, 

h2(nj,n2;m^,m2)  “  h(n  ^;m j)h(n 2:012)  (3.54) 

where  h(n;m)  =  C, ’sinclC^ ’(n-m)]  and  C,  and  are  constants.  Hence  »he  optical 
microscopic  imaging  system  of  the  square  incoherent  source  will  have  a  DIR  of  the  form 
in  (3.52). 


When  the  2-D  DIR  is  as  in  (3.52),  (3.51)  can  be  rewritten  as 
N-1  N-1  N-1  N-1 


V("i-n2)  -  2  ^_^Z^^^Z^^^q(nj;m^,ilj)q(n2,m2,l2)‘J(f^i-f"2)u(^i-^2) 

N-1  N-1  N-1  N-1  ,  , 

•  In?  ].„q(n2;m2,42)w<ni;m2l2) 


m2='0  2,2*0 


where 


j;m2,2,2)  A  2  ^  q(ni;m  1,2, )u(mi,m2)u(2 1,22)- 


(3.55) 


(3.56) 


The  minimum  norm  solution,  w(n j;m 2,2,2),  nj=0,1,*  • ‘.N-l,  can  be  obtained  via  the 
application  of  step  2)  in  the  1-D  algorithm.  Define, 


•  •  • 
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q(n  2,0,0)  q(n  2,0,1) 

q("2)  a  q(n2;1.0)  q(n2;l,l) 


q(n2;0,N-1) 

q(n,,1,N-l) 


02=0,1,*  •  •,N-1 


(3.57) 


|q(n2;N-l,0)q(n2;N-1,1)*  •  •  q(n2;N-1,N-l)J 

Then,  as  in  (3.45),  by  obtaining  the  dominant  eigensystem  of  the  Nx1  vectors  d„  , 

-"z 

n2*0,1,*  •  •,N-1,  are  obtained.  For  the  fixed  value  of  n^,  y(nj,n2)  can  be  obtained,  from 
dp^  and  w(n  ^;m2.1L2).  by 


''  N-1  N-1 

y(n2,n2)  =  ^  ^^d(n2;m2)d(n  2;!  2)^(0  i;m  2,4.2) 


(3.58) 


where  d(n;m)  is  the  m^^-element  value  of  the  vector  dp.  By  substituting  (3.56)  for 
w(n ^;m2.4.2),  (3.58)  can  be  rewritten  as,  for  n ^=0,1,*  •  •,N-1, 

«  N-1  N-1  N-1  N-1 

y(nj,n2)  •  t  A  .q{n  ,  ;m  ,  ,1 ,  )E  _l  _d(n  ,;m  ,)d(n  ,;8,  ,)u(m ,  ,m  ,)u(8, ,  ,8,  .13.59) 


N-1  N-1  N-1  N-1 


z(mi,n2)  4  ^  ^Qd(n2;m2)u(mj,m2).  (3.60) 

Then,  (3.59)  can  be  rewritten  as, 

A  N-1  N-1 

y(n^,n2)  •  Z  £  Q{n^;mi.li)z(m^.n^)z{l^,nA  (3.61) 

m  ^  W  Xr  ^  *w 

y(nj,n2)  is  given  in  (3.58)  and  q(n ^jm j,ll j)  is  also  known.  Hence,  (3.61)  is  the  1-0 
bilinear  expression  for  the  fixed  value  of  02,  02=0,1,*  •*, N-1.  By  applying  the  1-D 
algorithm  to  (3.61)  for  each  value  of  Oj,  n2=0.1,*  • ‘.N-l,  z(ftix,n2),  x  =0,1,*  •  •,N-1,  is 
obtained.  Then,  u(mj,m2)  can  be  solved  for  in  (3.60).  The  2-D  algorithm  with  the 
product  separabie  form  DIR  is  shown  on  Figure-7. 


The  algorithm  in  Figure-7  is  implemented  with 

h2(nj,n2;mj,m2)  »  sinc((ni-mJ/Nh]sincl(n2-m2)/Nh],  |nj-mj|<Nf,,  i»l,2  (3j62) 

0  otherwise 

where  sincx*[sin(trx)]/(TTx),  and 

y2(nii,m2, 11.1,4.2)  *  1  for  all  m  j,m 2.4 1,4. 2-  (3.63) 

Hence,  Figure-8(b)  shows  the  image  blurred  by  the  completely  coherent  optical  system 
with  Nji=4  from  the  31x31  original  image  in  Figure-8(a).  By  applying  the  algorithm  in 
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Figure-7,  the  image  shown  in  Figure-8(c)  has  been  obtained.  The  algorithm  results  in  the 
exact  image  restoration.  For  the  partially  coherent  image  restoration,  the  image  in 
Figure-9(a)  has  been  blurred  by  the  system  with  h2(n  ^,n  ^.m  ^,m  j)  in  (3.62)  and 

Yj{m^,m2,jli,i2)  »  sinc[{mj-aj)/Nj,]sinc((m2-2.2)/N^].  (3.64) 

The  blurred  image  is  shown  on  Figure-9(b).  has  been  chosen  as  19.  The  restored 
image  in  Figure-9(c)  retains  a  little  error  since  the  blurring  system  belongs  to  the  partially 
coherent  system  close  to  the  completely  coherent  case. 
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Obtain  Nxl  vectors  dp,  ,  n  * ‘.N-l.  by 

where  Xp  i  and  Vp  i  are  the  dominant 
eigenvaiue  and  eigenvector  of 


Obtain  the  minimum  norm  solution,  w(nj;m2,Z2),  of 
N-1  N-1 

y(n,,n2)  -  Z  Z  q(n2;m2.i2)w(np;m2,i2) 


m  2=012=0 


Compute 

V(ni<"2)  *  ^  lot 

where  d(n;m)  is  the  m^^-element  value  of  dp 


d{n  2;m  2)tl(n  ^ ;m  2,1 2) 


nj-n^+l 


n  j^N-l 


Yes  n2*0 


Obtain  2(mj,n2)  ^rom 

Y(n,.n2)  *  2  S  q(n  j;m  ^,1  ^)z(m  ^,n  2)z(2' i,n  2) 

^  *  m^»01^=0 

by  applying  step  2)  through  step  5)  of  1-0  algorithm. 


02=02+1 


n2=N-l 


I 


Yes  n^=0 


Solve  for  u(n  j,n2) 

2(01,02)  = 

N-1 

Z  d(n2;m2)u(ni.m2) 
m  2=0 

Figure-7.  2-0  Product  Separable  Form  Algorithm 
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4.  Image  Restoration  From  Noise  Added  Slurred  Image 


In  this  section  the  original  image  restoration  from  the  noise  added  image  following 
bilinear  blurring  in  (2.1)  with  nonnegative  OIR  will  be  considered.  The  blurred  image  with 
additive  noise  can  be  described  by 

g(n)  «  2  2  q(n;mj.m2)f(m.}f(mj)  +  v(n)  (4.1) 

where  v(n)  is  a  zero  mean  white  gaussian  noise.  The  bilinear  blurring  in  (4.1)  can  be 
recovered  by  applying  the  algorithm  in  Figure- 1.  But  due  to  the  noise  effect  the 
nonnegativity  condition  in  (2.4)  may  not  be  satisfied.  Then,  solving  (2.12)  recursively 
might  result  in  two  negative  solutions,  or  two  complex  solutions  in  addition  to  a 
nonnegative  solution  with  a  negative  one.  The  algorithm  in  Figure- 1  will  be  adapted  for 
use  in  this  case. 

Suppose  that  nonnegative  values  for  f(0),  f(1),*  *  *,  f(n-1)  have  been  obtained  and 

>  0,  i*1,2,***,n,  k=i,i+1,*  • ‘.N-l.  (4.2) 

Since  >  0  and  >  0  due  to  the  nonnegativity  of  q(n:mj,m2)  and 

Cn  -  Cf)  >  0  (4.3) 

due  to  (4.2),  solving  (2.12)  for  f(n)  will  result  in  one  nonnegative  solution,  say  s^,  and  one 
negative  solution,  Sj-  Let 

f^*Hn)  *  Sj.  (4.4) 

Using  f^^\n)  in  (4.4),  k»n+1,n+2.*  •  •,N-1  can  be  computed.  If 

>  0,  k-n+1,n+2,*  •  vN-l  (4.5) 

then 


f(n)  *  f(‘)(n)  (4.6) 

And  one  can  proceed  the  algorithm  to  solve  f(n+1). 

If  any  of  turns  out  to  be  negative,  say. 

ct"* ' )  <  0.  cp*  ‘ <  0,  •  •  •  ' )  <  0. 

then  by  choosing  as 


^  -f-  -- 


(4.7) 
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I 

[ 


I 


(4.8) 


the  output  index  I,  say  at  which  the  largest  modulus  occurs  among  negative 

•  •.fLj.  can  be  obtained.  For  to  be  nonnegative.  f{n)  has  to  satisfy 

-  [  q(il;n,n)f^(n)  +  2  ^  {q(2,;n,m)+q(il;m,n)}f(m)f(n)  ]  >  0  (4.9) 

*  m*0 

where  >  0  by  (4.3).  Let  0<a  <  1.  Then,  (4.9)  can  be  rewritten  by 

q(Jl;n,n)f ^(n)  +  Z  ^  {q(iL;n,m)+q(2,;m,n)}f(m)f{n)  -  aC^^»0.  (4.10) 

m»0  ** 

Choosing  proper  value  for  a  and  solving  (4.10)  for  f<n)  will  result  in  f^^Hn).  By  using  this 
the  nonnegativity  of  k»n+1,n+2.*  • ‘.N-l,  can  be  rechecked.  Note  that 


=  (1-a)c([’^  >  0. 


(4.11) 


If  alt  the  inequalities  in  (4.5)  are  satisfied  by  f^^^(n).  then  set 


f(n)  -  f<*)(n)  (4.12) 

2  %  •  and  proceed  the  algorithm  to  obtain  f(n-*-1).  If  any  one  or  more  of  inequalities  in  (4.5)  are 

not  satisfied  by  f^^Hn),  then,  by  repeating  the  above  procedure  until  all  the  inequalities  in 
(4.5)  can  be  satisfied,  one  can  obtain  a  desired  nonnegative  f(n).  It  is  easy  to  see  that  the 
number  of  negative  will  be  reduced  by,  at  least,  one  whenever  the  above 

I  procedure  is  repeated.  Hence,  by  repeating  the  above  procedure  at  most  r  times,  one  can 

obtain  a  nonnegative  f(n)  satisfying  (4.5).  This  scheme  is  shown  on  Figure- 10. 


In  the  simulation  SNR«30dB  white  gaussian  noise  is  added  to  the  blurred  image  in 
Figure-5(b).  The  resulting  noise  added  image  is  shown  on  Figure- 1 1(a).  By  choosing 
a«0.5,  the  restored  image  in  Figure-1 1(b)  is  obtained.  The  resulting  SNR  is  16dB. 


r 


r 
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Main  Algorithm 


Compute  ^ ) 

.  g(|f) 

(S|n-^0  -  c^n) 

- 

+  2  \q(k;m.n)-*-q{k;n,m)}f(m)f(n) 
m-0 

+  q(k;n.n)f^(n) 

n+1  <  k  <  N-1 

1 

Set 


cjn+0  >  0 


n+1  <  k  <  N-1 


1 
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10.  J.P.  Guiver  and  N.K.  Bose,  "Strictly  Hurwitz  property  invariance  of 
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14.  N.K.  Bose,  "Properties  of  Fade  approximants  to  Stieltjes  series  and 
systems,"  invited  paper  in  book  entitled,  "Rational  Approximation  and 
Interpolation" ,  ed.  by  P.R.  Graves-Morris ,  E.B.  Saff  and  R.S.  Varga, 
Springer- Verlag,  1985,  pp.  182-188. 

15.  N.K.  Bose,  "Applications  for  multidimensional  systems  theory,"  invited 
article  in  Encyclopedia  of  Systems  and  Control,  Pergamon  Press  Ltd., 
England,  scheduled  for  publication  in  1985. 

16.  N.K.  Bose,  "Multidimensional  Systems  stability,"  invited  article  in 
Encyclopedia  of  Systems  and  Control,  Pergamon  Press  Ltd.,  England, 
scheduled  for  publication  in  1985. 


17.  N\K.  Bose,  "Multivariate  realization  theory,"  invited  article  in  Ency¬ 
clopedia  of  Systems  and  Control,  Pergamon  Press  Ltd.,  England,  scheduled 
for  publication  in  1985. 
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e.  List  of  Personnel  Associated 


h  the  Research  Effort 


1.  N.K.  Bose 

Dr.  N.K.  Bose,  Professor  of  Electrical  Engineering  and  Professor  of  Mathe¬ 
matics  at  the  University  of  Pittsburgh  served  as  the  Principal  Investigator 
He  was  responsible  for  initiating  the  research,  conducting  it  to  its  success 
ful  completion  and  supervising  the  graduate  student . researchers  who  worked 
on  this  project. 

2.  Hector  M.  Valenzuela 

Served  as  a  graduate  student  researcher  for  eight  man-months,  from  March  1, 
1983  to  October  31,  1983.  He  completed  his  Ph.D.  dissertation  entitled, 
"Modeling  of  Multidimensional  Linear  Shift-Variant  Systems  with  Applications 
in  October  1983.  The  research  for  this  dissertation  was  supervised  by  Pro¬ 
fessor  N.K.  Bose. 

3 .  H . M .  Kim 

Served  as  a  graduate  student  researcher  for  13  man-months.  Mr.  Kim  is  cur¬ 
rently  a  candidate  for  the  Ph.D.  degree  in  Electrical  Engineering.  He  has 
been  performing  his  research  in  the  area  of  restoration  of  bilinearly  de¬ 
graded  images,  under  the  supervision  of  Professor  N.K.  Bose. 


f.  Interactions  of  Principal  Investigator 

.  Was  invited  to  be  a  member  of  the  Technical  Program  Committee  of  the  Inter¬ 
national  Symposium  on  Circuits  and  Systems  (ISCAS)  held  at  Newport  Beach, 
California,  May  2-4,  1983. 

Was  chairman  and  organizer  of  a  Special  Session  on  "Multidimensional  Systems: 
Spatio-Temporal  Filtering"  at  ISCAS,  Newport  Beach,  California,  May  2-4,  1983. 
).  Delivered  a  talk  on  "Modeling  of  2-D  LSV  systems  with  applications , "  based 
on  work  co-authored  with  Hector  M.  Valenzuela  at  ISCAS,  Newport  Beach,  Cali¬ 
fornia  May  2-4,  1983. 

Invited  to  give  a  seminar  on  "Rational  approximants  in  systems  theory,"  at 
the  U.K.  -  U,S.A.  Conference  on  Rational  Approximation  and  Interpolation, 
University  of  South  Florida,  Tampa,  Florida,  December  16,  1983. 

5.  Invited  to  give  a  seminar  on  "A  system-theoretic  approach  to  stability  of 
sets  of  polynomials,"  at  the  NSF  sponsored  Conference  on  "Linear  Algebra 
and  its  Role  in  Systems  Theory,"  Bowdoin  College,  Brunswick,  Maine, 

July  29- August  4,  1984. 

6.  Invited  to  give  a  talk  entitled,  "Novel  interpretations  in  multidimensional 
systems  theory,"  at  the  Special  Session  on  "Two-Dimensional  Circuits  and 
Systems  Theory  with  Applications,"  ISCAS,  Montreal,  May  10,  1984. 

7.  Invited  to  give  a  seminar  on  "Basic  tutorial  in  multidimensional  systems," 
at  Lehrstuhl  fur  Theoretische  Elektro technik  und  Messtechnik,  Universitaet 
Karlsruhe,  West  Germany,  July  19,  1984. 

8.  Invited  to  give  a  seminar  on  "Aspects  of  multidimensional  Systems  theory" 
at  Rutgers  University,  New  Brunswick,  New  Jersey,  on  March  21,  1985. 

9.  Invited  to  serve  as  a  member  of  program  committee  of  European  Conference 
on  Computer  Algebra  (Symbolic  and  Algebraic  Manipulation),  Linz,  Austria, 

April  1-3,  1985. 


10.  Invited  to  give  a  seminar  on  ’’Restoration  of  bilinearly  degraded  images, 
at  ISCAS,  Kyoto,  Japan,  June  7,  1985. 

11.  Invited  to  present  a  series  of  seminars  and  also  to  chair  a  session  on 
"Multidimensional  digital  signal  processing  2,"  at  the  Seventh  European 
Conference  on  Circuit  Theory  and  Design,  September  1985. 

12.  Invited  to  give  a  talk  on  "Status  of  recent  results  on  image  restoration 
in  Multidimensional  Systems  theory,"  at  the  24th  IEEE  Conference  on  Deci 
sion  and  Control,  Fort  Lauderdale,  Florida,  December  1985. 


jg»  Specific  Applications  Stemming  from  Research  Report 
Principal  Investigator:  Dr.  N.K.  Bose 

(a)  Primitive  Factorization  of  Bivariate  Polynomial  Matrices  Over  an  Arbitrary 
Field  of  Coefficients  I  Similar  to  the  widespread  use  of  irreducible  matrix 
fraction  descriptions  in  the  1-D  case,  such  representations  in  the  2-D  case 
are  known  to  have  g  "sat  potentials.  An  important  outcome  of  the  research 
conducted  has  been  the  presentation  of  a  primitive  factorization  theorem 
for  matrices  of  bivariate  polynomials  over  an  arbitrary  but  fixed  field  of 
coefficients  [1],  This,  in  turn,  leads  to  a  method  for  obtaining  an  ir¬ 
reducible  matrix  fraction  description  for  a  matrix  of  bivariate  rational 
functions  over  an  arbitrary  field.  In  fact  the  factorization  results  hold 
not  just  for  matrices  over  K[z,w]  (K  a  field),  but  for  matrices  over  D(w], 
where  D  is  an  arbitrary  Euclidean  domain  (or  in  theory,  over  a  Principal 
Ideal  Domain).  Importantly,  the  computations  required  to  obtain  the  fac¬ 
tors  in  K[z,w]  depend  neither  on  any  extension  field  nor  on  the  restric¬ 
tion  of  algebraic  closure,  in  contrast  to  earlier  approaches. 

The  matrix  fraction  descriptions  obtained  are  then  used  to  study 
stability  of  2-0  feedback  systems  where  the  plant  and  compensator  each 
corresponds  to  discrete  2-D  causal  or  weakly  causal  multi- input /multi¬ 
output  systems.  In  particular,  necessary  and  sufficient  conditions  are 
obtained  for  an  unstable  plant  to  be  stable  and  a  classification  of  the 
stabilizing  compensators  are  given  [2].  These  results  have  proven  and 
potential  applications  since  multidimensional  feedback  systems  have  been 
•  proposed  for  various  purposes  like  iterative  image  processing  and  restora¬ 
tion  13,4].  Such  image  processing  systems  that  contain  feedback  loops 
are  sometimes  known  to  oscillate  in  space  and  time  and  these  undesirable 
oscillations  can  only  be  avoided  if  proper  stability  conditions  are 
imposed  on  the  feedback  systems. 

(b)  Multidimensional  Linear  Shift-Variant  (LSV)  Systems:  For  any  n-D  discrete 
positive  cone  causal  (or  weakly  causal)  LSV  system,  whose  Impulse  response 
is  approximated  by  a  K-th  order  degenerate  sequence,  a  K-th  order  state- 
space  model  is  obtained  [5].  This  recursive  state-space  model  is  based  on 
a  (2”-l)-poincs  recurrence  formula,  which  for  the  causal  case  uses  the 
(2"-l)-cloaest neighboring  "past”  points  in  addition  to  the  input  in  order 

to  compute  any  current  output  state.  For  the  weakly  causal  case,  the  (2”-l) 
computed  points  required  are  not,  in  general,  the  closest  neighbors  to  the 
present  output,  which  is  being  computed.  Models  for  the  2-D  LSV  system  and 
its  inverse  can  be  used  to  perform  analysis  and  deconvolution  problems  very 
efficiently.  Examples  of  physically  motivated  applications  making  use  of 
the  theoretical  results  have  been  worked  out.  These  applications  include 
effects  of  1-0  LSV  motion  blur  and  the  blurring  due  to  Seidel  aberrations 
of  a  lens;  in  particular,  the  2-0  LSV  coma  aberration  was  studied  in  detail. 
The  reconstruction  of  the  original  object  from  the  LSV  blurred  image  was 
carried  out  successfully  by  means  of  the  state-space  model  for  the  inverse 
system.  There  are  several  advantages  of  the  approach  adopted  in  this  re¬ 
search.  First,  any  impulse  response  sequence  can  be  approximated  arbit¬ 
rarily  closely  by  a  K-th  order  degenerate  sequence  by  increasing  K.  The 
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recursive  iapiecencacion  is  associacec  vich  reduced  space- ric: 
conpucacionai  conplexity,  so  chac  ic  is  conceivable  that  degraded 
images  may  be  processed  in  real  time,  when  necessar/.  Though,  all 
types  of  point-spread  functions  cannot  be  modeled  via  the  recursive 
state-space  model,  the  flexibility  provided  by  the  weak  causality 
condition  broadens  considerably  the  scopes  for  applications.  Some¬ 
times  even  a  noncausal  point-spread  fimction  (as  illustrated  by  coma 
aberration  [5])  can  be  decomposed  suitably,  the  recursive  model 
applied  to  each  part,  and  the  results  carefully  superimposed  to 
yield  the  correct  solution.  Extensions  and  applicability  of  the 
results  to  image  restoration  subject  to  nonlinear  (especially 
bilinear)  degrading  phenomena  in  the  presence  of  nonstacionary 
noise  is  currently  being  investigated. 
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(c)  Restoration  of  Bilinearly  Degraded  Images:  Explained  in  adequate  detail  in 
Section  c  entitled  "Details  of  Research  Results  Obtained,"  of  this  report. 


